MAT8800 - Bayesian Regression Project
Duane Stanton
December 24, 2018
Introduction: Project, Resources, and Data
Project
The goal of this project is to apply the concepts of Bayesian regression as explored in my independent study on Bayesian regression in R and Stan, in coordination with Dr. Frey, during the Fall 2018 semester.
Given that the Fall term included consequential midterm Congressional elections (where party control of the House of Representatives and possibly the Senate were considered to possibly shift, and which did occur for the House), I chose to analyze House election data for the three most recent House elections - those for 2012, 2014, and 2016.
I determined House elections were a good candidate for such a regression analysis for several reasons:
The aforementioned overall topicality given we are in an “election year”.
With 435 Congressional districts per election, there is both a greater quantity of data and more granular data when considering the district level than comparative analysis for Senate elections (where senators serves for staggered 6-year terms such that 1/3 of Senate seats are “up for election” each election cycle - versus whole-House elections for 2-year terms for Representatives - and Senate elections occur at the state level, so there are fewer than 50 observations per election). More data and more granular data allows observed data to “have a greater say” in posterior distribution estimates relative to given prior estimates when fitting the models.
Because Congressional districts are required to be approximately equal in relative population size based on the most recently completed decennial Census, with a minimum of one Representative per state1, Congressional districts are by design less prone to having disproportionately large or small populations relative to other districts.
- Related to the previous point, the 2012, 2014, and 2016 elections each concerned districts based on population estimated from the 2010 Census, lending stability to district allocations.
Resources
This project was made possible thanks to several resources:
- Richard McElreath’s Statistical Rethinking2, the core text used for the independent study. It’s a great introduction to the overall mechanisms and basic theory of Bayesian regression, with most applications in base R and using the author’s accompanying
rethinkingR package.
- A Solomon Kurz’s Statistical Rethinking with brms, ggplot2, and the tidyverse3, a bookdown.org project working through the text using “tidyverse” rather than base R methodology, and the
brmspackage rather thanrethinking. This has been an excellent source of coding examples for working through McElreath’s text, with accompanying parallel insights and commentary.
- A. Vehtari, A. Gelman, and J. Gabry’s article Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC4, which provided insight into improving Bayesian model evaluation for out-of-sample prediction accuracy.
R tools including the R statistical programming language (v3.5.1 for the duration of my work), the R Studio integrated development environment (I used v1.1.456 for all my work, including creating this document), McElreath’s
rethinking(latest version 1.80), Paul Bürkner’sbrmspackage (Bayesian Regression Models using Stan, latest version 2.60), Hadley Wickham et. al.’s “tidyverse” suite of packages includingggplot2(v3.1.0) anddplyr(v0.7.8), and therstanpackage (v2.17.4), which interfaces with the Stan statistical modeling and computation platform which makes all the higher-level analysis possible.Additional R packages employed include
cowplot(v0.9.3) for side-by-side plots,bayesplot(v1.6.0) andtidybayes(v1.0.3) for Bayesian model diagnostics, andknitr(v1.2.0) to create this report as well as working withkableExtra(v0.9.0) to create summary tables in this report. I wholeheartedly apologize for any omissions of other packages used throughout this report, which should be displayed inlibrary(packageName)commands throughout.Last and certainly not least, Dr. Jesse Frey, who oversaw the independent study and provided greater insight into Bayesian methodology as well as advice on potential ways to strengthen this report and the overall independent study.
Please note: Any mistakes in this report are entirely my own!
Data
Data for the analysis comes from three sources:
- District-level election data comes from the Federal Election Commission’s (FEC’s) Election Results page5, for each year considered (2012, 2014, 2016). The data was read in from CSV files and processed to extract variables such as the per-district election winner and the maximum and total number of primary and general election candidates and votes by party (Democratic, Republican, and Other).
- Campaign-reported financial data comes from the FEC’s bulk downloads page of candidate financial summary data6 for each year considered - specifically the file
candidate_summary_[YYYY].csv. This data was matched to each campaign using the candidates’ FEC-assigned candidate ID, which was present in both the election and financial data sets.
Specifically, the variables for Total_Receipt (total campaign-reported money received, in USD) and Total_Disbursement (total campaign-reported money spent, also in USD) were used - this does not account for third-party spending (e.g. PACs and Super PACs), but provides some measure of reasonably well-defined financial data.
The top 10 entries based on Total_Disbursement for each year, and all entries with missing reported receipts and/or spending were investigated and confirmed against the FEC site and against the Center for Responsive Politics’ site7, which also monitors and records campaign-reported finances. This step turned up several cases of significant campaign activity (receipts/spending in excess of $100,000) as well as many smaller-level entries for each year that was not in the FEC bulk data.
- District-level demographics data comes from U.S. Census Bureau forms
B02001for ethnicity composition (counts of population identifying as White alone, Black or African American alone, American Indian and Alaska Native alone, Asian alone, Native Hawaiian and Other Pacific Islander alone, Some other race alone, or Two or more races), andS0201for all other demographics including male/female proportions, education level among the age 25+ population, and age-group proportions (among others). Data comes from American Community Survey (ACS) 1-year estimates for the year of each election considered, and were collected from the Census Bureau’s American FactFinder tool8.
Program Setup
# loading packages used
library(tidyverse)
library(knitr)
library(kableExtra)
library(brms)
library(bayesplot)
library(RColorBrewer)
# load 2012, 2014, 2016 Congressional District House elections and demographic data
houseElectionDat <-
read.csv("C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/houseElectionDat.csv")Considering the outcome variable - which party’s candidate was elected
While I also considered analyzing the vote margin between parties (e.g. # of votes for the Democratic candidate vs. # of votes for the Republican candidate, for each district), I settled on analyzing the outcome of a logistic regression model predicting whether or not a Democratic party candidate was elected for each district. This is primarily because a small set of districts had “uncontested” elections where either the Republican or Democratic candidate ran unopposed in the general election (and in a subset of these cases, was even unopposed in the primary).
At some future time, it may be interesting to consider a zero-inflated regression for the vote-margin model, however.
Winning party for each Congressional district
## # A tibble: 3 x 4
## year generalWinnerD generalWinnerR generalWinnerOth
## <int> <int> <int> <int>
## 1 2012 201 234 0
## 2 2014 188 247 0
## 3 2016 194 241 0
As the above plot show, Republican candidates won a majority of House seats in each election, with the greatest margin in 2014 (+59 seats). No third-party candidates won election in any race for the three years considered.
Because the Democratic party is the “challenging / minority” party for the 2018 midterm elections, I chose generalWinD (general election won by a Democratic candidate) as the outcome variable. Since no third party candidate was elected in the three years considered, the complement of this can be considered “general election won by a Republican candidate” for this analysis.
Considering predictor variables
Summary table of pre-standardization values for numeric predictor variables
The regression models below all use standardized continuous variables (with the exception of Boolean predictors, such as whether or not the [party] candidate is the incumbent, and small-count predictors, such as the number of primary candidates for [party]). This allows us to interpret posterior estimates for coefficients in terms of relative predictive association with a Democratic candidate being elected for relative changes in magnitude rather than in natural-scale units.
However, it is important to understand the natural scale of each variable, to put these magnitude changes in proper context.
| Predictor | Unit | Mean | SD | Mean | SD | Mean | SD |
|---|---|---|---|---|---|---|---|
| incumbentD | count | 3.890e-01 | 5.200e-01 | 4.210e-01 | 4.940e-01 | 3.950e-01 | 4.890e-01 |
| incumbentR | count | 5.150e-01 | 5.180e-01 | 4.800e-01 | 5.000e-01 | 4.920e-01 | 5.050e-01 |
| incumbentOth | count | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 |
| incumbentWinD | proportion | 3.490e-01 | 4.770e-01 | 3.930e-01 | 4.890e-01 | 3.840e-01 | 4.870e-01 |
| incumbentWinR | proportion | 4.570e-01 | 4.990e-01 | 4.670e-01 | 4.990e-01 | 4.710e-01 | 5.000e-01 |
| incumbentWinOth | proportion | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 |
| cmpgnTotRecptD | USD | 1.085e+06 | 1.273e+06 | 9.930e+05 | 1.231e+06 | 1.074e+06 | 1.642e+06 |
| cmpgnMaxRecptD | USD | 9.588e+05 | 1.044e+06 | 9.099e+05 | 1.069e+06 | 9.564e+05 | 1.293e+06 |
| cmpgnTotRecptR | USD | 1.391e+06 | 1.905e+06 | 1.325e+06 | 1.662e+06 | 1.260e+06 | 1.714e+06 |
| cmpgnMaxRecptR | USD | 1.250e+06 | 1.822e+06 | 1.171e+06 | 1.497e+06 | 1.120e+06 | 1.540e+06 |
| cmpgnTotRecptOth | USD | 3.039e+04 | 4.900e+05 | 8.299e+03 | 1.047e+05 | 2.522e+04 | 3.235e+05 |
| cmpgnMaxRecptOth | USD | 2.939e+04 | 4.726e+05 | 8.151e+03 | 1.035e+05 | 1.636e+04 | 1.791e+05 |
| cmpgnRecptTotal | USD | 2.507e+06 | 2.418e+06 | 2.326e+06 | 2.137e+06 | 2.360e+06 | 2.391e+06 |
| cmpgnTotDisbursD | USD | 1.073e+06 | 1.339e+06 | 9.381e+05 | 1.215e+06 | 9.602e+05 | 1.572e+06 |
| cmpgnMaxDisbursD | USD | 9.383e+05 | 1.047e+06 | 8.575e+05 | 1.052e+06 | 8.445e+05 | 1.213e+06 |
| cmpgnTotDisbursR | USD | 1.318e+06 | 1.828e+06 | 1.215e+06 | 1.628e+06 | 1.201e+06 | 1.582e+06 |
| cmpgnMaxDisbursR | USD | 1.178e+06 | 1.737e+06 | 1.064e+06 | 1.454e+06 | 1.062e+06 | 1.406e+06 |
| cmpgnMaxDisbursDvsR | USD | -2.394e+05 | 1.903e+06 | -2.067e+05 | 1.693e+06 | -2.179e+05 | 1.820e+06 |
| cmpgnTotDisbursOth | USD | 2.974e+04 | 4.881e+05 | 8.317e+03 | 1.046e+05 | 2.507e+04 | 3.227e+05 |
| cmpgnMaxDisbursOth | USD | 2.882e+04 | 4.720e+05 | 8.131e+03 | 1.034e+05 | 1.625e+04 | 1.784e+05 |
| cmpgnDisbursTotal | USD | 2.421e+06 | 2.423e+06 | 2.161e+06 | 2.154e+06 | 2.187e+06 | 2.337e+06 |
| primaryNumCandD | count | 1.924e+00 | 1.387e+00 | 1.621e+00 | 1.247e+00 | 1.775e+00 | 1.366e+00 |
| primaryTotD | vote count | 2.819e+04 | 2.601e+04 | 2.461e+04 | 2.351e+04 | 4.117e+04 | 4.096e+04 |
| primaryMaxD | vote count | 2.214e+04 | 2.170e+04 | 2.072e+04 | 2.069e+04 | 3.399e+04 | 3.596e+04 |
| primRunoffTotD | vote count | 4.704e+02 | 3.537e+03 | 9.917e+01 | 1.154e+03 | 9.303e+01 | 1.400e+03 |
| primRunoffMaxD | vote count | 2.794e+02 | 2.103e+03 | 5.660e+01 | 6.468e+02 | 5.541e+01 | 8.606e+02 |
| primaryUnopD | proportion | 1.680e-01 | 3.740e-01 | 1.790e-01 | 3.840e-01 | 1.860e-01 | 3.900e-01 |
| primaryNumCandR | count | 2.149e+00 | 1.753e+00 | 1.844e+00 | 1.516e+00 | 1.982e+00 | 1.833e+00 |
| primaryTotR | vote count | 3.613e+04 | 2.971e+04 | 2.937e+04 | 2.732e+04 | 4.404e+04 | 4.046e+04 |
| primaryMaxR | vote count | 2.740e+04 | 2.414e+04 | 2.200e+04 | 2.086e+04 | 3.416e+04 | 3.413e+04 |
| primRunoffTotR | vote count | 8.577e+02 | 5.720e+03 | 8.599e+02 | 6.348e+03 | 2.180e+02 | 3.039e+03 |
| primRunoffMaxR | vote count | 4.933e+02 | 3.255e+03 | 5.072e+02 | 3.803e+03 | 1.173e+02 | 1.635e+03 |
| primaryUnopR | proportion | 1.540e-01 | 3.610e-01 | 1.680e-01 | 3.740e-01 | 1.560e-01 | 3.640e-01 |
| primaryTotDvsR | vote count | -7.933e+03 | 3.932e+04 | -4.768e+03 | 3.604e+04 | -2.868e+03 | 5.517e+04 |
| primaryNumCandOth | count | 6.690e-01 | 1.016e+00 | 6.110e-01 | 8.950e-01 | 5.930e-01 | 8.970e-01 |
| primaryTotOth | vote count | 6.872e+02 | 2.912e+03 | 5.639e+02 | 2.786e+03 | 5.673e+02 | 2.634e+03 |
| primaryMaxOth | vote count | 6.138e+02 | 2.576e+03 | 4.809e+02 | 2.524e+03 | 4.537e+02 | 2.061e+03 |
| primRunoffTotOth | vote count | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 |
| primRunoffMaxOth | vote count | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 |
| primaryUnopOth | proportion | 4.160e-01 | 8.340e-01 | 4.070e-01 | 7.810e-01 | 4.620e-01 | 9.550e-01 |
| primaryNumCandTotal | count | 4.743e+00 | 2.637e+00 | 4.076e+00 | 2.356e+00 | 4.349e+00 | 2.782e+00 |
| primaryTotal | vote count | 6.501e+04 | 4.032e+04 | 5.454e+04 | 3.683e+04 | 8.578e+04 | 6.047e+04 |
| generalNumCandD | count | 9.890e-01 | 3.220e-01 | 9.910e-01 | 4.020e-01 | 1.039e+00 | 4.190e-01 |
| generalNumCandR | count | 1.014e+00 | 3.180e-01 | 9.750e-01 | 5.530e-01 | 9.930e-01 | 5.100e-01 |
| generalNumCandOth | count | 1.120e+00 | 1.071e+00 | 9.770e-01 | 1.058e+00 | 1.048e+00 | 1.334e+00 |
| generalUnopD | proportion | 2.000e-03 | 4.800e-02 | 2.000e-03 | 4.800e-02 | 2.000e-03 | 4.800e-02 |
| generalUnopR | proportion | 2.000e-03 | 4.800e-02 | 9.000e-03 | 9.600e-02 | 2.000e-03 | 4.800e-02 |
| generalUnopOth | proportion | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 |
| generalWinD | proportion | 4.620e-01 | 4.990e-01 | 4.320e-01 | 4.960e-01 | 4.460e-01 | 4.980e-01 |
| generalWinR | proportion | 5.380e-01 | 4.990e-01 | 5.680e-01 | 4.960e-01 | 5.540e-01 | 4.980e-01 |
| generalWinOth | proportion | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 | 0.000e+00 |
| generalNumCandTotal | count | 3.122e+00 | 1.128e+00 | 2.943e+00 | 1.246e+00 | 3.080e+00 | 1.486e+00 |
| generalTotal | vote count | 2.789e+05 | 6.152e+04 | 1.792e+05 | 5.832e+04 | 2.959e+05 | 6.213e+04 |
| totalPop | count | 7.202e+05 | 3.482e+04 | 7.315e+05 | 3.882e+04 | 7.413e+05 | 4.631e+04 |
| popDensityPerSqMi | count per sq. mi. | 2.369e+03 | 6.885e+03 | 2.418e+03 | 7.051e+03 | 2.440e+03 | 7.115e+03 |
| pctMale | percentage | 4.920e+01 | 9.610e-01 | 4.921e+01 | 9.640e-01 | 4.923e+01 | 9.520e-01 |
| pctFemale | percentage | 5.080e+01 | 9.610e-01 | 5.079e+01 | 9.640e-01 | 5.077e+01 | 9.520e-01 |
| medianAgeYr | years | 3.767e+01 | 3.647e+00 | 3.799e+01 | 3.659e+00 | 3.829e+01 | 3.668e+00 |
| age18_34Pop | count | 1.686e+05 | 2.558e+04 | 1.716e+05 | 2.630e+04 | 1.728e+05 | 2.714e+04 |
| age35_64Pop | count | 2.834e+05 | 2.134e+04 | 2.849e+05 | 2.129e+04 | 2.865e+05 | 2.288e+04 |
| age65plusPop | count | 9.901e+04 | 2.282e+04 | 1.061e+05 | 2.369e+04 | 1.130e+05 | 2.516e+04 |
| pctEduHSorLessAge25plus | percentage | 4.191e+01 | 9.511e+00 | 4.110e+01 | 9.497e+00 | 4.002e+01 | 9.397e+00 |
| pctEduSomColgAge25plus | percentage | 2.931e+01 | 4.486e+00 | 2.920e+01 | 4.552e+00 | 2.909e+01 | 4.653e+00 |
| pctEduBachAge25plus | percentage | 1.804e+01 | 5.670e+00 | 1.847e+01 | 5.725e+00 | 1.913e+01 | 5.832e+00 |
| pctEduGradProAge25plus | percentage | 1.074e+01 | 4.783e+00 | 1.123e+01 | 4.956e+00 | 1.176e+01 | 5.101e+00 |
| pctEduBachPlusAge25plus | percentage | 2.877e+01 | 1.008e+01 | 2.969e+01 | 1.032e+01 | 3.089e+01 | 1.055e+01 |
| pctUSBorn | percentage | 8.703e+01 | 1.107e+01 | 8.681e+01 | 1.104e+01 | 8.660e+01 | 1.119e+01 |
| pctForgnBornUSCtzn | percentage | 5.937e+00 | 5.634e+00 | 1.319e+01 | 1.104e+01 | 1.340e+01 | 1.119e+01 |
| pctAge5PlusNonEnglshHome | percentage | 2.109e+01 | 1.780e+01 | 2.107e+01 | 1.779e+01 | 2.148e+01 | 1.781e+01 |
| pctInLaborForceUnemp | percentage | 5.953e+00 | 1.648e+00 | 4.555e+00 | 1.259e+00 | 3.616e+00 | 9.370e-01 |
| pctNotInLaborForce | percentage | 3.617e+01 | 4.809e+00 | 3.674e+01 | 4.996e+00 | 3.691e+01 | 5.082e+00 |
| medianHHIncome | USD | 5.330e+04 | 1.427e+04 | 5.595e+04 | 1.513e+04 | 6.005e+04 | 1.638e+04 |
| perCapitaIncome | USD | 2.726e+04 | 7.520e+03 | 2.881e+04 | 8.097e+03 | 3.104e+04 | 8.637e+03 |
| pctNoHealthInsCivNonInst | percentage | 1.478e+01 | 5.787e+00 | 1.166e+01 | 4.872e+00 | 8.532e+00 | 4.330e+00 |
| pctPovertyAge18_64 | percentage | 1.485e+01 | 5.135e+00 | 1.464e+01 | 5.056e+00 | 1.327e+01 | 4.544e+00 |
| pctWhiteAlone | percentage | 7.405e+01 | 1.711e+01 | 7.356e+01 | 1.717e+01 | 7.281e+01 | 1.724e+01 |
| pctAfrAmAlone | percentage | 1.250e+01 | 1.427e+01 | 1.255e+01 | 1.421e+01 | 1.256e+01 | 1.391e+01 |
| pctOtherAlone | percentage | 1.056e+01 | 1.032e+01 | 1.089e+01 | 1.064e+01 | 1.141e+01 | 1.111e+01 |
| pctMultiRacial | percentage | 2.888e+00 | 1.903e+00 | 3.005e+00 | 1.964e+00 | 3.215e+00 | 2.008e+00 |
| log_popDensPerSqMi | log(count per sq. mi.) | 6.043e+00 | 1.896e+00 | 6.059e+00 | 1.899e+00 | 6.071e+00 | 1.900e+00 |
Evaluating possible predictor variables - EDA and thoughts on potential posteriors
Given the abundance of predictors which could be included in the model, I decided to focus on those which I felt were most likely to have some meaningful (very plausibly non-zero) association with the predictor variable.
Additionally, I wanted predictor variables to provide good breadth of coverage for what I feel are certain key factors relevant to understanding election outcomes:
Incumbency - whether a given district has the previously-elected candidate in the running, and - if so - which party they are a member of
Campaign finances - some measure of the resources available (or utilized) for election contestants
Primary election results - some measure of political activity in the primary elections to choose general-election candidates
Demographic factors - things such as population density (to get an idea of relative urbanicity/rurality), age or educational mix, and economic status, among others
Summary of district-level predictors considered
incumbentDandincumbentR- district has an incumbent for D/R party contesting the election (count)cmpgnMaxDisbursDvsR- difference in maximum reported campaign spend by D/R candidate (I considered using TOTAL spend by party here, but correlation between max and total spend is very high here - \(\approx 0.95\) for D and \(\approx 0.97\) for R across the full dataset)primaryTotDvsR- difference in total vote count for each party in primary elections (may consider primary winner’s share of total party primary votes in future iteration)primaryUnopDandprimaryUnopR- district has a D/R candidate running unopposed in the primary (Boolean)popDensityPerSqMi- district population density per square milemedianAgeYr- district median agepctEduHSorLessAge25plus- district % of age 25+ population with HS or less as highest completed education levelpctEduGradProAge25plus- district % of age 25+ population with Graduate/Professional degree as highest completed education levelpctInLaborForceUnemp- % of district labor force (age 16+, either in civilian labor force (employed or seeking work, not in an institution e.g. school or prison) OR armed forces) unemployedperCapitaIncome- total district income divided by total district populationpctWhiteAlone- % of total district population identifying as ‘white’ only
Potential alternative/additional predictors OR for another model - for future consideration:
cmpgnTotDisbursD and cmpgnTotDisbursR, primaryMax[D/R] / primaryTot[D/R], generalNumCandTotal, pctFemale, age18_34Pop and age65plusPop, pctEduGradProAge25plus, pctUSBorn AND/OR pctForgnBornUSCtzn, pctNotInLaborForce, pctNoHealthInsCivNonInst, pctPovertyAge18_64, pctAfrAmAlone, pctMultiRacial
Exploratory Data Analysis - predictors and outcome variables
While these plots use non-standardized data for the sake of exploring the data in each variable’s natural scale, the fit models employ standardized numerical variables to allow greater ease of interpreting predictor coefficients in terms of relative magnitudes.
Note: Each predictor will be assigned a moderately regularizing normal prior, centered at 0.
Incumbency and outcome
The above shows that, in the vast majority of cases where there is an incumbent (standing for re-election) from a given party, that incumbent wins. I anticipate a relatively strongly positive posterior distribution for incumbentD and an equally negative posterior for incumbentR.
Reported campaign spending - maximum vs. total by party (per district)
Maximum and Total spend by party / year / district are generally quite closely related - for now I’ll consider maximum campaign spend levels.
Reported campaign spending - maximum D minus R vs. outcome by party (per district)
It appears that districts where one party has greater reported spend relative to the other is more likely to win that election. It seems reasonable to expect a weakly positive posterior in this case.
Primary votes - maximum vs. total by party (per district)
While this is somewhat similar to the comparison of maximum / total reported campaign spend, the greater fan-shaped distribution of points indicates there is somewhat less of a correlation between maximum / total primary votes (by party) than in the case of maximum / total campaign spend.
I may be introducing inconsistency in my approach here, but I feel that maximum campaign spend is a more relevant metric than total spend, while the reverse is true in terms of primary votes. This is because, I believe, that total primary vote count for a given party is indicative of the enthusiasm/interest for a given party in a district. While maximum primary vote count may reflect the winning candidate’s strength relative to their competition, I anticipate that voters are more likely to back their preferred party and to some degree disregard the primary winner’s relative strength than they are to put greater weight on primary win margins.
As with difference in maximum campaign spend, difference in total primary vote by party seems to provide a reasonably good indicator of the ultimate likelihood a Democratic candidate wins election in the House.
It seems prudent to compare the two predictors against each other, in case there is excessive correlation between the two variables:
While there is certainly evidence of a positive correlation between the party difference in maximum campaign spend and the party difference in total primary votes, it’s not worryingly high:
houseElectionDat %>%
group_by(year) %>%
summarize(corVal = round(cor((cmpgnMaxDisbursD - cmpgnMaxDisbursR),
(primaryTotD - primaryTotR)), 3) ) %>%
kable() %>%
kable_styling(full_width = F)| year | corVal |
|---|---|
| 2012 | 0.461 |
| 2014 | 0.463 |
| 2016 | 0.530 |
I’m not sure, but the widening discrepancy over time may be reflective of increasing partisan lean within districts over time. It should be noted that the 2016 presidential campaign was cast as especially partisan, however, so it may be related to that.
Primary unopposed D/R candidate
The plot above suggests a candidates’ running unopposed in their party’s primary election doesn’t seems to provide a great deal of insight into the likelihood they (or at least their party) win the general election. Also, in the majority of elections, both the Republican and Democratic primaries have multiple challengers.
Population density
Most districts have relatively low density, but Democratic candidates fare disproportionately well in high-density districts. A moderately weakly positive posterior seems a reasonable expectation.
Because log-population density has a much more symmetrical distribution, I’ll use that. [Updated 12/22/18]
Median age
It seems that Democratic candidates might fare relatively better in districts with lower median ages, and Republican candidates fare better in districts with higher median ages, but there’s nothing especially impressive here.
Perhaps looking at district composition by proportions in younger / older age groups would be useful?
As with the median age plot, there’s mildly suggestive evidence here, but nothing conclusive. I’ll stick with medianAgeYr for now, and expect a weak, mildly negative posterior.
Education level among the age 25+ population
I originally also considered using “undergraduate or greater” education here, but there was no meaningful pattern there, so I revised to “graduate/professional degree”.
There is marginal at best evidence of HS or less education having greater association with Democratic or Republican electoral success, but there is reasonably good evidence of districts with greater concentration of graduate/professional degree holders favoring Democratic candidates.
I’ll drop the HS or less predictor, and anticipate a weakly positive posterior for the graduate/professional degree predictor.
Labor force unemployment rate
There is some evidence of districts with higher unemployment rates having a greater association with Democratic candidate election winners. Here it seems that a weakly positive posterior could arise.
Let’s check to what extent unemployment rate appears to be correlated with graduate/professional education:
While the districts with the greatest relative proportions of graduate/professional degree-holders have consistently low unemployment rates, there is increasingly great variability in unemployment rates as the district proportion of graduate/professional degree-holders decreases. This results in a weak negative correlation between the two variables.
Per capita income
There also appears to be a positive association between greater per capita income for a district and that district’s electing a Democratic candidate. Again, a weakly positive posterior seems plausible.
As with unemployment rate, it seems prudent to consider the extent to which per capita income is correlated with the proportion of highly credentialed inhabitants, and also with the unemployment rate.
There is a strong, consistent positive correlation between graduate/professional degree and per capita income. If the inclusion of both predictors becomes problematic for fitting the model, I’m inclined to drop the per capita income predictor, as this would seem to be the more difficult of the two to accurately measure in a “real world” setting faced by a non-governmental entity interested in estimating House election outcomes.
There is also a weaker, consistent negative correlation between unemployment rate and per capita income.
Proportion of population identifying as white
Districts where 50% or less of the population identifies as white only appear to only elect Democratic candidates in the three years considered (with the exception of a single district in 2016). For districts where more than half of the population identifies as white only, Democratic candidates are elected, but not nearly as often as are Republican candidates.
I anticipate a negative posterior will emerge here.
Model fitting and analysis
My original model-fitting approach was the following:
Fit models for standardized 2012 data only, using varying subsets of predictors, and compare models by information criterion values and by information criterion-based relative model weights. Also confirm the models were fit well using trace plots and related diagnostics.
Among the “top-performing” model(s) selected, evaluate in-sample prediction performance, and consider whether trimming the “full” model (with all considered predictors) leads to improved IC metrics.
For the same model(s) chosen from step 2, make predictions using 2014 data and evaluate.
Fit the “full” model to 2014 standardized data and assess model fit as well as whether a “reduced” version has superior IC metrics, and any major differences compared to the 2012 “full” model. Repeat assessment of performance in predicting “in sample” (2014) and then “out of sample” (2016) data.
Repeat step 4 for 2016 data, without the “out of sample” prediction.
Fit “all years full” models, this time considering fixed-effects (intercept and slopes) vs. varying-effects models since now there are 3 observations for each district (one for each year’s election). Analyze the fit and in-sample performance of these models and conclude that the two models are not meaningfully different, and relegate the analysis to the Appendix.
Fit “all years full” models with binomial-outcome variables (aggregating predictors across the three years considered for each district). Analyze the fit and in-sample performance of these models and go forward with sensitivity analysis (considering the effect of looser priors on posterior estimates) and counterfactual analysis (holding other predictors constant at defined levels, assessing the impact on predicted outcomes when manipulating the level of a single predictor) for the “top” model.
**In order to “get to the good stuff in this report, step 7 from above is shown in the next section, and all work from steps 1 through 5 is in Appendix 1, while all work from step 6 is in Appendix 2.**
Fitting and Analyzing “all years” models
# load standardized predictors set
houseElectionDat_std <-
read.csv("C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/houseElectionDat_std.csv")
# break out 2012, 2014, 2016 data
houseElectionDat12_std <-
houseElectionDat_std %>%
filter(year == 2012)
houseElectionDat14_std <-
houseElectionDat_std %>%
filter(year == 2014)
houseElectionDat16_std <-
houseElectionDat_std %>%
filter(year == 2016)Model 2 - Binomial model, using 2012 / 2014 / 2016 data
This section treats the outcome of a given districts set of elections as a Binomial random variable (with the outcome being a Democratic candidate’s election 0, 1, 2, or 3 times out of 3 elections).
Rather than treating each year’s district-level election as unique events, let’s fit fixed-effects and varying-effects models evaluating each district’s propensity to elect a Democratic or Republican candidate - in other words, aggregating data at the district level. This should also hopefully provide a proper setting for partial pooling to improve the model fit.
To aggregate predictor data, I’ll take the mean of the three years’ observations for each predictor. Note: While I also tried fitting a model incorporating measurement error (as the standard error of each district’s mean predictor value), it is apparently more computationally demanding than I’d thought, so I’m sticking with the two models shown below.
# first need to aggregate predictor estimates and response outcomes at district level
# here I take the mean of each district-level predictor
# and include district-level measurement error as sd(predictor) / sqrt(nObs = 3)
houseElectionDat_std_agg <-
houseElectionDat_std %>%
group_by(district2) %>%
# only going to retain predictors previously used in full models thus far
# note: aside from the outcome var, predictors are MEAN values
# note: adding a very small positive constant to _SE vars to resolve issue
summarize(nObs = n(),
generalWinD = sum(generalWinD),
incumbentD_mn = mean(incumbentD),
incumbentD_SE = sd(incumbentD) / sqrt(nObs),
incumbentR_mn = mean(incumbentR),
incumbentR_SE = sd(incumbentR) / sqrt(nObs),
primaryUnopD_mn = mean(primaryUnopD),
primaryUnopD_SE = sd(primaryUnopD) / sqrt(nObs),
primaryUnopR_mn = mean(primaryUnopR),
primaryUnopR_SE = sd(primaryUnopR) / sqrt(nObs),
cmpgnMaxDisbursDvsR_mn = mean(cmpgnMaxDisbursDvsR),
cmpgnMaxDisbursDvsR_SE = sd(cmpgnMaxDisbursDvsR) / sqrt(nObs),
primaryTotDvsR_mn = mean(primaryTotDvsR),
primaryTotDvsR_SE = sd(primaryTotDvsR) / sqrt(nObs),
popDensityPerSqMi_mn = mean(popDensityPerSqMi),
popDensityPerSqMi_SE = sd(popDensityPerSqMi) / sqrt(nObs),
medianAgeYr_mn = mean(medianAgeYr),
medianAgeYr_SE = sd(medianAgeYr) / sqrt(nObs),
pctEduGradProAge25plus_mn = mean(pctEduGradProAge25plus),
pctEduGradProAge25plus_SE = sd(pctEduGradProAge25plus) / sqrt(nObs),
pctInLaborForceUnemp_mn = mean(pctInLaborForceUnemp),
pctInLaborForceUnemp_SE = sd(pctInLaborForceUnemp) / sqrt(nObs),
perCapitaIncome_mn = mean(perCapitaIncome),
perCapitaIncome_SE = sd(perCapitaIncome) / sqrt(nObs),
pctWhiteAlone_mn = mean(pctWhiteAlone),
pctWhiteAlone_SE = sd(pctWhiteAlone) / sqrt(nObs),
log_popDensPerSqMi_mn = mean(log_popDensPerSqMi),
log_popDensPerSqMi_SE = sd(log_popDensPerSqMi) / sqrt(nObs) )Fitting the models, and preliminary fit summaries
m2.allFull_fIfS <-
brm(generalWinD | trials(nObs) ~ 1 + incumbentD_mn + incumbentR_mn +
primaryUnopD_mn + primaryUnopR_mn +
cmpgnMaxDisbursDvsR_mn + primaryTotDvsR_mn +
log_popDensPerSqMi_mn + medianAgeYr_mn +
pctEduGradProAge25plus_mn +
pctInLaborForceUnemp_mn + perCapitaIncome_mn +
pctWhiteAlone_mn,
data = houseElectionDat_std_agg, family = binomial(),
prior = c(set_prior("normal(0, 3)", class = "Intercept"),
set_prior("normal(0, 2)", class = "b")),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99),
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m2.allFull_fIfS")
m2.allFull_vIvS <-
brm(generalWinD | trials(nObs) ~ 1 + incumbentD_mn + incumbentR_mn +
primaryUnopD_mn + primaryUnopR_mn + cmpgnMaxDisbursDvsR_mn +
primaryTotDvsR_mn + log_popDensPerSqMi_mn + medianAgeYr_mn +
pctEduGradProAge25plus_mn + pctInLaborForceUnemp_mn +
perCapitaIncome_mn + pctWhiteAlone_mn +
(1 + incumbentD_mn + incumbentR_mn +
primaryUnopD_mn + primaryUnopR_mn + cmpgnMaxDisbursDvsR_mn +
primaryTotDvsR_mn + log_popDensPerSqMi_mn + medianAgeYr_mn +
pctEduGradProAge25plus_mn + pctInLaborForceUnemp_mn +
perCapitaIncome_mn + pctWhiteAlone_mn |district2),
data = houseElectionDat_std_agg, family = binomial(),
prior = c(set_prior("normal(0, 3)", class = "Intercept"),
set_prior("normal(0, 2)", class = "b"),
set_prior("cauchy(0, 2)", class = "sd"),
set_prior("lkj(2)", class = "cor")),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99),
thin = 2,
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m2.allFull_vIvS")summary(m2.allFull_fIfS)## Family: binomial
## Links: mu = logit
## Formula: generalWinD | trials(nObs) ~ 1 + incumbentD_mn + incumbentR_mn + primaryUnopD_mn + primaryUnopR_mn + cmpgnMaxDisbursDvsR_mn + primaryTotDvsR_mn + log_popDensPerSqMi_mn + medianAgeYr_mn + pctEduGradProAge25plus_mn + pctInLaborForceUnemp_mn + perCapitaIncome_mn + pctWhiteAlone_mn
## Data: houseElectionDat_std_agg (Number of observations: 435)
## Samples: 3 chains, each with iter = 7000; warmup = 3000; thin = 1;
## total post-warmup samples = 12000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept -0.59 0.75 -2.07 0.86 12000 1.00
## incumbentD_mn 4.76 0.91 3.01 6.60 10733 1.00
## incumbentR_mn -2.93 0.97 -4.85 -1.05 10777 1.00
## primaryUnopD_mn 0.27 0.84 -1.40 1.90 9076 1.00
## primaryUnopR_mn -0.51 0.90 -2.26 1.28 8968 1.00
## cmpgnMaxDisbursDvsR_mn 1.82 0.46 0.95 2.76 12000 1.00
## primaryTotDvsR_mn 1.06 0.50 0.09 2.05 12000 1.00
## log_popDensPerSqMi_mn 0.28 0.31 -0.31 0.88 12000 1.00
## medianAgeYr_mn -0.11 0.26 -0.61 0.39 12000 1.00
## pctEduGradProAge25plus_mn -0.48 0.56 -1.59 0.63 9099 1.00
## pctInLaborForceUnemp_mn 0.46 0.41 -0.36 1.26 12000 1.00
## perCapitaIncome_mn 1.08 0.63 -0.16 2.31 7776 1.00
## pctWhiteAlone_mn -0.60 0.45 -1.52 0.24 9908 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).The fixed-effect model has good effective sample size (maximum possible: n = 12,000) and convergence diagnostics (R-hat = exactly 1). There is certainly variety in the predictor posterior distribution estimates, with some having 95% credible intervals not including 0 and others essentially centered over 0.
summary(m2.allFull_vIvS)## Family: binomial
## Links: mu = logit
## Formula: generalWinD | trials(nObs) ~ 1 + incumbentD_mn + incumbentR_mn + primaryUnopD_mn + primaryUnopR_mn + cmpgnMaxDisbursDvsR_mn + primaryTotDvsR_mn + log_popDensPerSqMi_mn + medianAgeYr_mn + pctEduGradProAge25plus_mn + pctInLaborForceUnemp_mn + perCapitaIncome_mn + pctWhiteAlone_mn + (1 + incumbentD_mn + incumbentR_mn + primaryUnopD_mn + primaryUnopR_mn + cmpgnMaxDisbursDvsR_mn + primaryTotDvsR_mn + log_popDensPerSqMi_mn + medianAgeYr_mn + pctEduGradProAge25plus_mn + pctInLaborForceUnemp_mn + perCapitaIncome_mn + pctWhiteAlone_mn | district2)
## Data: houseElectionDat_std_agg (Number of observations: 435)
## Samples: 3 chains, each with iter = 7000; warmup = 3000; thin = 2;
## total post-warmup samples = 6000
##
## Group-Level Effects:
## ~district2 (Number of levels: 435)
## Estimate Est.Error l-95% CI
## sd(Intercept) 0.47 0.36 0.02
## sd(incumbentD_mn) 0.84 0.64 0.03
## sd(incumbentR_mn) 0.89 0.64 0.04
## sd(primaryUnopD_mn) 1.33 0.99 0.05
## sd(primaryUnopR_mn) 0.90 0.73 0.04
## sd(cmpgnMaxDisbursDvsR_mn) 0.69 0.53 0.03
## sd(primaryTotDvsR_mn) 0.75 0.57 0.03
## sd(log_popDensPerSqMi_mn) 0.46 0.36 0.02
## sd(medianAgeYr_mn) 0.69 0.42 0.04
## sd(pctEduGradProAge25plus_mn) 0.43 0.33 0.02
## sd(pctInLaborForceUnemp_mn) 0.50 0.39 0.02
## sd(perCapitaIncome_mn) 0.49 0.37 0.02
## sd(pctWhiteAlone_mn) 0.64 0.47 0.03
## cor(Intercept,incumbentD_mn) -0.01 0.25 -0.49
## cor(Intercept,incumbentR_mn) -0.02 0.25 -0.50
## cor(incumbentD_mn,incumbentR_mn) 0.00 0.25 -0.47
## cor(Intercept,primaryUnopD_mn) -0.00 0.25 -0.47
## cor(incumbentD_mn,primaryUnopD_mn) 0.00 0.25 -0.48
## cor(incumbentR_mn,primaryUnopD_mn) 0.01 0.25 -0.47
## cor(Intercept,primaryUnopR_mn) -0.01 0.25 -0.49
## cor(incumbentD_mn,primaryUnopR_mn) -0.01 0.25 -0.49
## cor(incumbentR_mn,primaryUnopR_mn) -0.01 0.25 -0.48
## cor(primaryUnopD_mn,primaryUnopR_mn) -0.02 0.25 -0.52
## cor(Intercept,cmpgnMaxDisbursDvsR_mn) -0.01 0.25 -0.49
## cor(incumbentD_mn,cmpgnMaxDisbursDvsR_mn) -0.02 0.25 -0.50
## cor(incumbentR_mn,cmpgnMaxDisbursDvsR_mn) 0.00 0.25 -0.47
## cor(primaryUnopD_mn,cmpgnMaxDisbursDvsR_mn) -0.02 0.25 -0.50
## cor(primaryUnopR_mn,cmpgnMaxDisbursDvsR_mn) -0.00 0.25 -0.47
## cor(Intercept,primaryTotDvsR_mn) -0.01 0.25 -0.48
## cor(incumbentD_mn,primaryTotDvsR_mn) -0.02 0.25 -0.50
## cor(incumbentR_mn,primaryTotDvsR_mn) -0.00 0.25 -0.49
## cor(primaryUnopD_mn,primaryTotDvsR_mn) -0.02 0.25 -0.50
## cor(primaryUnopR_mn,primaryTotDvsR_mn) -0.00 0.25 -0.50
## cor(cmpgnMaxDisbursDvsR_mn,primaryTotDvsR_mn) -0.02 0.25 -0.50
## cor(Intercept,log_popDensPerSqMi_mn) 0.02 0.25 -0.47
## cor(incumbentD_mn,log_popDensPerSqMi_mn) 0.00 0.25 -0.48
## cor(incumbentR_mn,log_popDensPerSqMi_mn) 0.03 0.26 -0.46
## cor(primaryUnopD_mn,log_popDensPerSqMi_mn) 0.03 0.25 -0.46
## cor(primaryUnopR_mn,log_popDensPerSqMi_mn) -0.01 0.25 -0.48
## cor(cmpgnMaxDisbursDvsR_mn,log_popDensPerSqMi_mn) -0.01 0.25 -0.50
## cor(primaryTotDvsR_mn,log_popDensPerSqMi_mn) -0.03 0.25 -0.50
## cor(Intercept,medianAgeYr_mn) 0.01 0.25 -0.48
## cor(incumbentD_mn,medianAgeYr_mn) 0.02 0.25 -0.46
## cor(incumbentR_mn,medianAgeYr_mn) 0.01 0.24 -0.46
## cor(primaryUnopD_mn,medianAgeYr_mn) 0.01 0.25 -0.47
## cor(primaryUnopR_mn,medianAgeYr_mn) -0.00 0.25 -0.47
## cor(cmpgnMaxDisbursDvsR_mn,medianAgeYr_mn) 0.01 0.25 -0.47
## cor(primaryTotDvsR_mn,medianAgeYr_mn) -0.00 0.25 -0.47
## cor(log_popDensPerSqMi_mn,medianAgeYr_mn) 0.01 0.25 -0.47
## cor(Intercept,pctEduGradProAge25plus_mn) 0.00 0.25 -0.48
## cor(incumbentD_mn,pctEduGradProAge25plus_mn) 0.00 0.25 -0.47
## cor(incumbentR_mn,pctEduGradProAge25plus_mn) -0.00 0.25 -0.47
## cor(primaryUnopD_mn,pctEduGradProAge25plus_mn) 0.01 0.25 -0.48
## cor(primaryUnopR_mn,pctEduGradProAge25plus_mn) -0.00 0.25 -0.49
## cor(cmpgnMaxDisbursDvsR_mn,pctEduGradProAge25plus_mn) 0.00 0.25 -0.47
## cor(primaryTotDvsR_mn,pctEduGradProAge25plus_mn) -0.00 0.25 -0.48
## cor(log_popDensPerSqMi_mn,pctEduGradProAge25plus_mn) -0.02 0.25 -0.51
## cor(medianAgeYr_mn,pctEduGradProAge25plus_mn) 0.00 0.25 -0.48
## cor(Intercept,pctInLaborForceUnemp_mn) -0.01 0.24 -0.47
## cor(incumbentD_mn,pctInLaborForceUnemp_mn) -0.02 0.25 -0.49
## cor(incumbentR_mn,pctInLaborForceUnemp_mn) 0.00 0.25 -0.48
## cor(primaryUnopD_mn,pctInLaborForceUnemp_mn) -0.01 0.25 -0.50
## cor(primaryUnopR_mn,pctInLaborForceUnemp_mn) -0.00 0.25 -0.48
## cor(cmpgnMaxDisbursDvsR_mn,pctInLaborForceUnemp_mn) -0.00 0.25 -0.49
## cor(primaryTotDvsR_mn,pctInLaborForceUnemp_mn) -0.00 0.25 -0.47
## cor(log_popDensPerSqMi_mn,pctInLaborForceUnemp_mn) -0.01 0.25 -0.49
## cor(medianAgeYr_mn,pctInLaborForceUnemp_mn) -0.00 0.25 -0.49
## cor(pctEduGradProAge25plus_mn,pctInLaborForceUnemp_mn) 0.02 0.25 -0.47
## cor(Intercept,perCapitaIncome_mn) -0.00 0.25 -0.50
## cor(incumbentD_mn,perCapitaIncome_mn) 0.01 0.25 -0.48
## cor(incumbentR_mn,perCapitaIncome_mn) -0.00 0.25 -0.49
## cor(primaryUnopD_mn,perCapitaIncome_mn) 0.02 0.25 -0.47
## cor(primaryUnopR_mn,perCapitaIncome_mn) 0.00 0.25 -0.47
## cor(cmpgnMaxDisbursDvsR_mn,perCapitaIncome_mn) 0.01 0.25 -0.47
## cor(primaryTotDvsR_mn,perCapitaIncome_mn) 0.00 0.25 -0.48
## cor(log_popDensPerSqMi_mn,perCapitaIncome_mn) -0.01 0.25 -0.49
## cor(medianAgeYr_mn,perCapitaIncome_mn) 0.00 0.25 -0.47
## cor(pctEduGradProAge25plus_mn,perCapitaIncome_mn) -0.03 0.26 -0.53
## cor(pctInLaborForceUnemp_mn,perCapitaIncome_mn) 0.01 0.25 -0.47
## cor(Intercept,pctWhiteAlone_mn) 0.00 0.25 -0.48
## cor(incumbentD_mn,pctWhiteAlone_mn) 0.01 0.25 -0.47
## cor(incumbentR_mn,pctWhiteAlone_mn) 0.01 0.25 -0.48
## cor(primaryUnopD_mn,pctWhiteAlone_mn) 0.02 0.25 -0.45
## cor(primaryUnopR_mn,pctWhiteAlone_mn) 0.01 0.25 -0.47
## cor(cmpgnMaxDisbursDvsR_mn,pctWhiteAlone_mn) 0.00 0.25 -0.48
## cor(primaryTotDvsR_mn,pctWhiteAlone_mn) 0.01 0.25 -0.48
## cor(log_popDensPerSqMi_mn,pctWhiteAlone_mn) 0.04 0.25 -0.47
## cor(medianAgeYr_mn,pctWhiteAlone_mn) 0.00 0.25 -0.47
## cor(pctEduGradProAge25plus_mn,pctWhiteAlone_mn) 0.01 0.25 -0.48
## cor(pctInLaborForceUnemp_mn,pctWhiteAlone_mn) 0.01 0.25 -0.47
## cor(perCapitaIncome_mn,pctWhiteAlone_mn) 0.01 0.25 -0.47
## u-95% CI Eff.Sample Rhat
## sd(Intercept) 1.33 4286 1.00
## sd(incumbentD_mn) 2.40 5080 1.00
## sd(incumbentR_mn) 2.35 3842 1.00
## sd(primaryUnopD_mn) 3.70 3970 1.00
## sd(primaryUnopR_mn) 2.71 4988 1.00
## sd(cmpgnMaxDisbursDvsR_mn) 1.97 4952 1.00
## sd(primaryTotDvsR_mn) 2.11 4642 1.00
## sd(log_popDensPerSqMi_mn) 1.35 4978 1.00
## sd(medianAgeYr_mn) 1.61 3517 1.00
## sd(pctEduGradProAge25plus_mn) 1.25 5295 1.00
## sd(pctInLaborForceUnemp_mn) 1.43 5294 1.00
## sd(perCapitaIncome_mn) 1.38 5176 1.00
## sd(pctWhiteAlone_mn) 1.72 4171 1.00
## cor(Intercept,incumbentD_mn) 0.48 4390 1.00
## cor(Intercept,incumbentR_mn) 0.46 4667 1.00
## cor(incumbentD_mn,incumbentR_mn) 0.48 4918 1.00
## cor(Intercept,primaryUnopD_mn) 0.48 4553 1.00
## cor(incumbentD_mn,primaryUnopD_mn) 0.48 4640 1.00
## cor(incumbentR_mn,primaryUnopD_mn) 0.48 4584 1.00
## cor(Intercept,primaryUnopR_mn) 0.47 4344 1.00
## cor(incumbentD_mn,primaryUnopR_mn) 0.48 4205 1.00
## cor(incumbentR_mn,primaryUnopR_mn) 0.47 4446 1.00
## cor(primaryUnopD_mn,primaryUnopR_mn) 0.45 4368 1.00
## cor(Intercept,cmpgnMaxDisbursDvsR_mn) 0.48 4275 1.00
## cor(incumbentD_mn,cmpgnMaxDisbursDvsR_mn) 0.46 4495 1.00
## cor(incumbentR_mn,cmpgnMaxDisbursDvsR_mn) 0.49 4443 1.00
## cor(primaryUnopD_mn,cmpgnMaxDisbursDvsR_mn) 0.46 4558 1.00
## cor(primaryUnopR_mn,cmpgnMaxDisbursDvsR_mn) 0.48 4255 1.00
## cor(Intercept,primaryTotDvsR_mn) 0.47 4491 1.00
## cor(incumbentD_mn,primaryTotDvsR_mn) 0.46 4531 1.00
## cor(incumbentR_mn,primaryTotDvsR_mn) 0.47 4517 1.00
## cor(primaryUnopD_mn,primaryTotDvsR_mn) 0.47 4715 1.00
## cor(primaryUnopR_mn,primaryTotDvsR_mn) 0.49 4752 1.00
## cor(cmpgnMaxDisbursDvsR_mn,primaryTotDvsR_mn) 0.46 4461 1.00
## cor(Intercept,log_popDensPerSqMi_mn) 0.50 4615 1.00
## cor(incumbentD_mn,log_popDensPerSqMi_mn) 0.49 4470 1.00
## cor(incumbentR_mn,log_popDensPerSqMi_mn) 0.52 4894 1.00
## cor(primaryUnopD_mn,log_popDensPerSqMi_mn) 0.49 4861 1.00
## cor(primaryUnopR_mn,log_popDensPerSqMi_mn) 0.48 4632 1.00
## cor(cmpgnMaxDisbursDvsR_mn,log_popDensPerSqMi_mn) 0.47 3970 1.00
## cor(primaryTotDvsR_mn,log_popDensPerSqMi_mn) 0.46 4244 1.00
## cor(Intercept,medianAgeYr_mn) 0.49 5027 1.00
## cor(incumbentD_mn,medianAgeYr_mn) 0.50 4938 1.00
## cor(incumbentR_mn,medianAgeYr_mn) 0.47 4455 1.00
## cor(primaryUnopD_mn,medianAgeYr_mn) 0.50 4951 1.00
## cor(primaryUnopR_mn,medianAgeYr_mn) 0.47 4909 1.00
## cor(cmpgnMaxDisbursDvsR_mn,medianAgeYr_mn) 0.50 4710 1.00
## cor(primaryTotDvsR_mn,medianAgeYr_mn) 0.48 5190 1.00
## cor(log_popDensPerSqMi_mn,medianAgeYr_mn) 0.49 4817 1.00
## cor(Intercept,pctEduGradProAge25plus_mn) 0.49 4435 1.00
## cor(incumbentD_mn,pctEduGradProAge25plus_mn) 0.48 4411 1.00
## cor(incumbentR_mn,pctEduGradProAge25plus_mn) 0.47 4180 1.00
## cor(primaryUnopD_mn,pctEduGradProAge25plus_mn) 0.49 3936 1.00
## cor(primaryUnopR_mn,pctEduGradProAge25plus_mn) 0.48 4351 1.00
## cor(cmpgnMaxDisbursDvsR_mn,pctEduGradProAge25plus_mn) 0.47 4342 1.00
## cor(primaryTotDvsR_mn,pctEduGradProAge25plus_mn) 0.48 4181 1.00
## cor(log_popDensPerSqMi_mn,pctEduGradProAge25plus_mn) 0.46 3976 1.00
## cor(medianAgeYr_mn,pctEduGradProAge25plus_mn) 0.48 4755 1.00
## cor(Intercept,pctInLaborForceUnemp_mn) 0.46 4338 1.00
## cor(incumbentD_mn,pctInLaborForceUnemp_mn) 0.48 4545 1.00
## cor(incumbentR_mn,pctInLaborForceUnemp_mn) 0.49 4416 1.00
## cor(primaryUnopD_mn,pctInLaborForceUnemp_mn) 0.48 4530 1.00
## cor(primaryUnopR_mn,pctInLaborForceUnemp_mn) 0.48 4857 1.00
## cor(cmpgnMaxDisbursDvsR_mn,pctInLaborForceUnemp_mn) 0.49 4390 1.00
## cor(primaryTotDvsR_mn,pctInLaborForceUnemp_mn) 0.47 4223 1.00
## cor(log_popDensPerSqMi_mn,pctInLaborForceUnemp_mn) 0.48 4677 1.00
## cor(medianAgeYr_mn,pctInLaborForceUnemp_mn) 0.47 4891 1.00
## cor(pctEduGradProAge25plus_mn,pctInLaborForceUnemp_mn) 0.50 4493 1.00
## cor(Intercept,perCapitaIncome_mn) 0.48 4441 1.00
## cor(incumbentD_mn,perCapitaIncome_mn) 0.49 4390 1.00
## cor(incumbentR_mn,perCapitaIncome_mn) 0.49 4811 1.00
## cor(primaryUnopD_mn,perCapitaIncome_mn) 0.50 4395 1.00
## cor(primaryUnopR_mn,perCapitaIncome_mn) 0.48 4351 1.00
## cor(cmpgnMaxDisbursDvsR_mn,perCapitaIncome_mn) 0.49 4425 1.00
## cor(primaryTotDvsR_mn,perCapitaIncome_mn) 0.49 3268 1.00
## cor(log_popDensPerSqMi_mn,perCapitaIncome_mn) 0.47 4683 1.00
## cor(medianAgeYr_mn,perCapitaIncome_mn) 0.49 4745 1.00
## cor(pctEduGradProAge25plus_mn,perCapitaIncome_mn) 0.46 4343 1.00
## cor(pctInLaborForceUnemp_mn,perCapitaIncome_mn) 0.49 4510 1.00
## cor(Intercept,pctWhiteAlone_mn) 0.49 4381 1.00
## cor(incumbentD_mn,pctWhiteAlone_mn) 0.49 4369 1.00
## cor(incumbentR_mn,pctWhiteAlone_mn) 0.47 4683 1.00
## cor(primaryUnopD_mn,pctWhiteAlone_mn) 0.50 4587 1.00
## cor(primaryUnopR_mn,pctWhiteAlone_mn) 0.50 4590 1.00
## cor(cmpgnMaxDisbursDvsR_mn,pctWhiteAlone_mn) 0.49 4403 1.00
## cor(primaryTotDvsR_mn,pctWhiteAlone_mn) 0.49 4530 1.00
## cor(log_popDensPerSqMi_mn,pctWhiteAlone_mn) 0.51 3794 1.00
## cor(medianAgeYr_mn,pctWhiteAlone_mn) 0.49 5061 1.00
## cor(pctEduGradProAge25plus_mn,pctWhiteAlone_mn) 0.48 4890 1.00
## cor(pctInLaborForceUnemp_mn,pctWhiteAlone_mn) 0.49 4473 1.00
## cor(perCapitaIncome_mn,pctWhiteAlone_mn) 0.49 4724 1.00
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept -0.40 1.00 -2.34 1.61 5058 1.00
## incumbentD_mn 5.63 1.25 3.24 8.13 5248 1.00
## incumbentR_mn -3.75 1.31 -6.37 -1.25 4936 1.00
## primaryUnopD_mn 0.48 1.11 -1.69 2.63 5316 1.00
## primaryUnopR_mn -1.12 1.16 -3.44 1.08 5298 1.00
## cmpgnMaxDisbursDvsR_mn 2.92 0.79 1.49 4.57 4840 1.00
## primaryTotDvsR_mn 1.85 0.78 0.40 3.47 5119 1.00
## log_popDensPerSqMi_mn 0.56 0.47 -0.31 1.52 5228 1.00
## medianAgeYr_mn -0.13 0.38 -0.89 0.60 5152 1.00
## pctEduGradProAge25plus_mn -0.06 0.79 -1.59 1.51 5414 1.00
## pctInLaborForceUnemp_mn 0.60 0.59 -0.56 1.73 4986 1.00
## perCapitaIncome_mn 0.76 0.82 -0.86 2.37 5634 1.00
## pctWhiteAlone_mn -0.82 0.64 -2.12 0.39 4842 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).The varying-effects model is decidedly more complex, with sd parameters estimating the variability of partially pooled predictor posteriors (these seem to indicate there is appreciable variability for each), cor parameters estimating the correlation of posterior slopes (I think that, due to the brms default of using non-centered parametrization, it’s unsurprising these are all very plausibly 0), and population-level posteriors which appear to be relatively different from that of the fixed-effects model.
To get a better idea of this last set, let’s compare the population-level predictor posteriors between the fixed- and varying-effects models:
While the varying-effects model has wider 90% credible intervals (indicating greater uncertainty in the posterior estimates) for each population-level parameter, it does not have all point estimates shrinking towards zero as I might have expected. Indeed, for parameters whose fixed-effects estimates are further from zero, the corresponding varying-effects estimates are further still from zero - though the two 90% credible intervals largely overlap in all cases.
Model comparisons - relative IC values/weights, and in-sample prediction performance
Relative Information Criterion values
I anticipate that the greater ability of the varying-effects model to “learn from the data” (in the sense of potentially partially pool intercept/slope estimates across multiple, similar districts) will result in relatively superior IC values and, correspondingly, a greater share of IC-allocated model weight.
(m2.allFull_kFold <- kfold(m2.allFull_fIfS, m2.allFull_vIvS,
K = 10, compare = T) )## KFOLDIC SE
## m2.allFull_fIfS 160.60 27.42
## m2.allFull_vIvS 160.17 25.02
## m2.allFull_fIfS - m2.allFull_vIvS 0.43 7.87(m2.allFull_modWts <-
model_weights(m2.allFull_fIfS, m2.allFull_vIvS, weights = "kfold", K = 10) )## m2.allFull_fIfS m2.allFull_vIvS
## 0.7006 0.2994K-fold cross-validation IC-based weighting assigns a roughly 70/30 split between the fixed-effects and varying-effects models.
In-sample prediction performance
Let’s assess how well m2.allFull_fIfS and m2.allFull_vIvS predict the number of Democratic candidates elected across the three elections, and how the two models compare (for in-sample data).
# prediction of each district (columns) for each post-warmup iteration (rows)
postPred_m2.allFull_fIfS <- posterior_predict(m2.allFull_fIfS)
# mean for each per-year district (column of postPred_m2.allFull_fIfS)
meanPred_m2.allFull_fIfS <- apply(postPred_m2.allFull_fIfS, 2, mean)
predVsObs.2f <-
tibble(district2 = houseElectionDat_std_agg$district2,
generalWinD = houseElectionDat_std_agg$generalWinD,
meanPredWinD = meanPred_m2.allFull_fIfS )
nBins <- 100
cols <- c("darkred", "red",
"lightgrey",
"blue", "darkblue")
colGradient <- colorRampPalette(cols)
cut.cols <- colGradient(50)
cuts <- cut(predVsObs.2f$meanPredWinD, nBins)
names(cuts) <- sapply(cuts,function(t)
cut.cols[which(as.character(t) == levels(cuts))])
predVsObs.2f %>%
ggplot(aes(x = meanPredWinD, fill = cut(meanPredWinD, nBins))) +
geom_histogram(bins = nBins) +
facet_wrap(. ~ generalWinD) +
labs(title = "Posterior prediction check for m2.allFull_fIfS",
subtitle = "Faceted by # Democratic candidates elected (out of 3 elections: 2012/14/16)",
x = "Mean count of model predictions, D candidate elected",
y = "Count") +
scale_fill_manual(values = cut.cols,
guide = F) +
theme_ds1()
# prediction of each district (columns) for each post-warmup iteration (rows)
postPred_m2.allFull_vIvS <- posterior_predict(m2.allFull_vIvS)
# mean for each per-year district (column of postPred_m2.allFull_vIvS)
meanPred_m2.allFull_vIvS <- apply(postPred_m2.allFull_vIvS, 2, mean)
predVsObs.2v <-
tibble(district2 = houseElectionDat_std_agg$district2,
generalWinD = houseElectionDat_std_agg$generalWinD,
meanPredWinD = meanPred_m2.allFull_vIvS )
nBins <- 100
cols <- c("darkred", "red",
"lightgrey",
"blue", "darkblue")
colGradient <- colorRampPalette(cols)
cut.cols <- colGradient(50)
cuts <- cut(predVsObs.2v$meanPredWinD, nBins)
names(cuts) <- sapply(cuts,function(t)
cut.cols[which(as.character(t) == levels(cuts))])
predVsObs.2v %>%
ggplot(aes(x = meanPredWinD, fill = cut(meanPredWinD, nBins))) +
geom_histogram(bins = nBins) +
facet_wrap(. ~ generalWinD) +
labs(title = "Posterior prediction check for m2.allFull_vIvS",
subtitle = "Faceted by # Democratic candidates elected (out of 3 elections: 2012/14/16)",
x = "Mean count of model predictions, D candidate elected",
y = "Count") +
scale_fill_manual(values = cut.cols,
guide = F) +
theme_ds1()
# comparison of fixed- and varying-effects model predictions
ggplot() +
geom_density(data = predVsObs.2f,
aes(x = meanPredWinD), fill = color_set8[5], alpha = 0.7,
adjust = 0.1) +
geom_density(data = predVsObs.2v,
aes(x = meanPredWinD), fill = color_set8[7], alpha = 0.7,
adjust = 0.1) +
labs(title = "Comparing posterior prediction densities, 2016 models",
subtitle = "Fixed (tan) vs Varying effects (orange)",
x = "Mean proportion of model predictions, D candidate elected",
y = "Density") +
theme_ds1()
While there are relatively few observations with Democratic candidates being elected in 1 or 2 out of the three elections (in most cases the districts are either “fully Republican” or “fully Democrat”), the varying-effects model seems to outperform the fixed-effects model in these cases, as well as in the “all or nothing” outcomes which make up the majority.
In order to compare the two models, let’s define a “miss” as a model predicting a mean outcome (# of Democratic candidates elected among the 3 elections for each district) that is 0.5 or more away from the observed outcome.
predVsObs.2 <-
tibble(district2 = houseElectionDat_std_agg$district2,
generalWinD = houseElectionDat_std_agg$generalWinD,
meanPredWinD.f = meanPred_m2.allFull_fIfS,
meanPredWinD.v = meanPred_m2.allFull_vIvS) %>%
mutate(hitOrMiss.f = if_else((generalWinD == 0 & meanPredWinD.f < 0.5) |
(generalWinD == 1 & meanPredWinD.f > 0.5 &
meanPredWinD.f < 1.5) |
(generalWinD == 2 & meanPredWinD.f > 1.5 &
meanPredWinD.f < 2.5) |
(generalWinD == 3 & meanPredWinD.f > 2.5 &
meanPredWinD.f < 3.5) |
(generalWinD == 4 & meanPredWinD.f > 3.5),
"hit", "miss"),
hitOrMiss.v = if_else((generalWinD == 0 & meanPredWinD.v < 0.5) |
(generalWinD == 1 & meanPredWinD.v > 0.5 &
meanPredWinD.v < 1.5) |
(generalWinD == 2 & meanPredWinD.v > 1.5 &
meanPredWinD.v < 2.5) |
(generalWinD == 3 & meanPredWinD.v > 2.5 &
meanPredWinD.v < 3.5) |
(generalWinD == 4 & meanPredWinD.v > 3.5),
"hit", "miss") )
predVsObs.2 %>%
group_by(generalWinD) %>%
summarize(hit.f = sum(hitOrMiss.f == "hit"),
miss.f = sum(hitOrMiss.f == "miss"),
hit.v = sum(hitOrMiss.v == "hit"),
miss.v = sum(hitOrMiss.v == "miss") ) %>%
kable(align = "c") %>%
kable_styling(full_width = F, bootstrap_options = "striped")| generalWinD | hit.f | miss.f | hit.v | miss.v |
|---|---|---|---|---|
| 0 | 222 | 3 | 225 | 0 |
| 1 | 8 | 13 | 13 | 8 |
| 2 | 1 | 4 | 4 | 1 |
| 3 | 183 | 1 | 183 | 1 |
As expected from the above plots, the varying-effects model “misses” less frequently for each observed outcome count. The fixed-effects model “misses” in 21 cases, while the varying-effects model “misses” in 11 cases - and for each level of generalWinD, the varying-effects model has at most half as many “misses” as seen in the fixed-effects case.
This is counter to what I would have expected based on K-fold IC weights.
For either model, the “misses” are disproportionately more likely in cases where a Democratic candidate was actually elected once out of the three possible elections. For three of the four potential outcomes (where a Democratic candidate is elected 0, 1, or 2 times across the three elections), the varying-effects model is less likely to “miss” than the fixed-effects model.
Therefore, I will go forward with the varying-effects model for further analysis, including considering possible interactions.
Sensitivity analysis
As discussed with Dr. Frey, here I consider the extent to which relaxing prior distribution estimates impacts the posterior distribution of the model. The following pair of models use weaker priors which impart less regularization of parameter estimates compared to the prior pair of fixed- and varying-effects models.
m2.allFull_fIfS2 <-
brm(generalWinD | trials(nObs) ~ 1 + incumbentD_mn + incumbentR_mn +
primaryUnopD_mn + primaryUnopR_mn +
cmpgnMaxDisbursDvsR_mn + primaryTotDvsR_mn +
log_popDensPerSqMi_mn + medianAgeYr_mn +
pctEduGradProAge25plus_mn +
pctInLaborForceUnemp_mn + perCapitaIncome_mn +
pctWhiteAlone_mn,
data = houseElectionDat_std_agg, family = binomial(),
prior = c(set_prior("normal(0, 10)", class = "Intercept"),
set_prior("normal(0, 10)", class = "b")),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99),
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m2.allFull_fIfS2")
m2.allFull_vIvS2 <-
brm(generalWinD | trials(nObs) ~ 1 + incumbentD_mn + incumbentR_mn +
primaryUnopD_mn + primaryUnopR_mn + cmpgnMaxDisbursDvsR_mn +
primaryTotDvsR_mn + log_popDensPerSqMi_mn + medianAgeYr_mn +
pctEduGradProAge25plus_mn + pctInLaborForceUnemp_mn +
perCapitaIncome_mn + pctWhiteAlone_mn +
(1 + incumbentD_mn + incumbentR_mn +
primaryUnopD_mn + primaryUnopR_mn + cmpgnMaxDisbursDvsR_mn +
primaryTotDvsR_mn + log_popDensPerSqMi_mn + medianAgeYr_mn +
pctEduGradProAge25plus_mn + pctInLaborForceUnemp_mn +
perCapitaIncome_mn + pctWhiteAlone_mn |district2),
data = houseElectionDat_std_agg, family = binomial(),
prior = c(set_prior("normal(0, 10)", class = "Intercept"),
set_prior("normal(0, 10)", class = "b"),
set_prior("cauchy(0, 10)", class = "sd"),
set_prior("lkj(1)", class = "cor")),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99),
thin = 2,
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m2.allFull_vIvS2")Let’s check the model summaries to review posterior estimates and the R-hat values:
summary(m2.allFull_fIfS2)## Family: binomial
## Links: mu = logit
## Formula: generalWinD | trials(nObs) ~ 1 + incumbentD_mn + incumbentR_mn + primaryUnopD_mn + primaryUnopR_mn + cmpgnMaxDisbursDvsR_mn + primaryTotDvsR_mn + log_popDensPerSqMi_mn + medianAgeYr_mn + pctEduGradProAge25plus_mn + pctInLaborForceUnemp_mn + perCapitaIncome_mn + pctWhiteAlone_mn
## Data: houseElectionDat_std_agg (Number of observations: 435)
## Samples: 3 chains, each with iter = 7000; warmup = 3000; thin = 1;
## total post-warmup samples = 12000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept -0.90 0.96 -2.81 0.95 10006 1.00
## incumbentD_mn 5.75 1.20 3.45 8.18 9706 1.00
## incumbentR_mn -3.00 1.24 -5.53 -0.60 9376 1.00
## primaryUnopD_mn 0.27 1.06 -1.86 2.28 8407 1.00
## primaryUnopR_mn -0.34 1.16 -2.65 1.89 7896 1.00
## cmpgnMaxDisbursDvsR_mn 1.76 0.51 0.82 2.80 12000 1.00
## primaryTotDvsR_mn 0.98 0.57 -0.10 2.10 9764 1.00
## log_popDensPerSqMi_mn 0.25 0.33 -0.38 0.91 12000 1.00
## medianAgeYr_mn -0.15 0.28 -0.69 0.38 10940 1.00
## pctEduGradProAge25plus_mn -0.80 0.62 -2.04 0.44 9496 1.00
## pctInLaborForceUnemp_mn 0.54 0.45 -0.35 1.40 10173 1.00
## perCapitaIncome_mn 1.51 0.72 0.13 2.94 8178 1.00
## pctWhiteAlone_mn -0.67 0.49 -1.68 0.23 9187 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).The posterior estimates seem reasonable (no wildly extended ranges), effective sample sizes are fairly large, and the R-hat values suggest convergence has been achieved. Overlaid density plots and trace plots also confirmed good model gits (posterior estimate density plots aligned across each MCMC chain’s estimates, and each chain’s trace plot line indicated good stationarity over a relatively narrow range of values, and good mixing of chains).
Let’s compare posterior distribution estimates for the “more regularized” versus “less regularized” fixed-effects models:
The only parameter which appears to be much different between the fixed-effects models is that the posterior for incumbentD appears greater for the “less regularized” model (though a majority of the 90% credible intervals overlap). Credible intervals are also a bit wider for the “less regularized” model.
Now let’s consider the “less regularized” varying-effects model:
summary(m2.allFull_vIvS2)## Family: binomial
## Links: mu = logit
## Formula: generalWinD | trials(nObs) ~ 1 + incumbentD_mn + incumbentR_mn + primaryUnopD_mn + primaryUnopR_mn + cmpgnMaxDisbursDvsR_mn + primaryTotDvsR_mn + log_popDensPerSqMi_mn + medianAgeYr_mn + pctEduGradProAge25plus_mn + pctInLaborForceUnemp_mn + perCapitaIncome_mn + pctWhiteAlone_mn + (1 + incumbentD_mn + incumbentR_mn + primaryUnopD_mn + primaryUnopR_mn + cmpgnMaxDisbursDvsR_mn + primaryTotDvsR_mn + log_popDensPerSqMi_mn + medianAgeYr_mn + pctEduGradProAge25plus_mn + pctInLaborForceUnemp_mn + perCapitaIncome_mn + pctWhiteAlone_mn | district2)
## Data: houseElectionDat_std_agg (Number of observations: 435)
## Samples: 3 chains, each with iter = 7000; warmup = 3000; thin = 2;
## total post-warmup samples = 6000
##
## Group-Level Effects:
## ~district2 (Number of levels: 435)
## Estimate Est.Error l-95% CI
## sd(Intercept) 1.32 1.04 0.05
## sd(incumbentD_mn) 2.56 1.91 0.11
## sd(incumbentR_mn) 3.12 1.99 0.17
## sd(primaryUnopD_mn) 3.91 2.72 0.17
## sd(primaryUnopR_mn) 2.52 2.01 0.10
## sd(cmpgnMaxDisbursDvsR_mn) 1.72 1.41 0.07
## sd(primaryTotDvsR_mn) 2.22 1.69 0.08
## sd(log_popDensPerSqMi_mn) 1.32 1.03 0.05
## sd(medianAgeYr_mn) 1.46 1.05 0.07
## sd(pctEduGradProAge25plus_mn) 1.06 0.87 0.04
## sd(pctInLaborForceUnemp_mn) 1.29 1.07 0.04
## sd(perCapitaIncome_mn) 1.26 0.98 0.05
## sd(pctWhiteAlone_mn) 1.75 1.29 0.07
## cor(Intercept,incumbentD_mn) -0.02 0.27 -0.53
## cor(Intercept,incumbentR_mn) -0.04 0.27 -0.56
## cor(incumbentD_mn,incumbentR_mn) -0.01 0.26 -0.52
## cor(Intercept,primaryUnopD_mn) 0.01 0.26 -0.49
## cor(incumbentD_mn,primaryUnopD_mn) -0.01 0.27 -0.53
## cor(incumbentR_mn,primaryUnopD_mn) 0.03 0.27 -0.50
## cor(Intercept,primaryUnopR_mn) -0.02 0.27 -0.52
## cor(incumbentD_mn,primaryUnopR_mn) -0.01 0.27 -0.52
## cor(incumbentR_mn,primaryUnopR_mn) -0.03 0.27 -0.54
## cor(primaryUnopD_mn,primaryUnopR_mn) -0.05 0.27 -0.57
## cor(Intercept,cmpgnMaxDisbursDvsR_mn) -0.01 0.27 -0.51
## cor(incumbentD_mn,cmpgnMaxDisbursDvsR_mn) -0.03 0.27 -0.53
## cor(incumbentR_mn,cmpgnMaxDisbursDvsR_mn) 0.01 0.27 -0.50
## cor(primaryUnopD_mn,cmpgnMaxDisbursDvsR_mn) -0.03 0.26 -0.53
## cor(primaryUnopR_mn,cmpgnMaxDisbursDvsR_mn) -0.00 0.27 -0.51
## cor(Intercept,primaryTotDvsR_mn) -0.01 0.27 -0.52
## cor(incumbentD_mn,primaryTotDvsR_mn) -0.01 0.26 -0.52
## cor(incumbentR_mn,primaryTotDvsR_mn) -0.00 0.26 -0.50
## cor(primaryUnopD_mn,primaryTotDvsR_mn) -0.03 0.26 -0.54
## cor(primaryUnopR_mn,primaryTotDvsR_mn) -0.00 0.27 -0.51
## cor(cmpgnMaxDisbursDvsR_mn,primaryTotDvsR_mn) -0.02 0.27 -0.54
## cor(Intercept,log_popDensPerSqMi_mn) 0.03 0.27 -0.50
## cor(incumbentD_mn,log_popDensPerSqMi_mn) 0.01 0.27 -0.50
## cor(incumbentR_mn,log_popDensPerSqMi_mn) 0.05 0.27 -0.48
## cor(primaryUnopD_mn,log_popDensPerSqMi_mn) 0.05 0.27 -0.48
## cor(primaryUnopR_mn,log_popDensPerSqMi_mn) -0.01 0.26 -0.51
## cor(cmpgnMaxDisbursDvsR_mn,log_popDensPerSqMi_mn) -0.01 0.27 -0.52
## cor(primaryTotDvsR_mn,log_popDensPerSqMi_mn) -0.05 0.27 -0.54
## cor(Intercept,medianAgeYr_mn) 0.00 0.26 -0.51
## cor(incumbentD_mn,medianAgeYr_mn) 0.01 0.27 -0.51
## cor(incumbentR_mn,medianAgeYr_mn) 0.00 0.26 -0.49
## cor(primaryUnopD_mn,medianAgeYr_mn) -0.01 0.26 -0.51
## cor(primaryUnopR_mn,medianAgeYr_mn) -0.01 0.26 -0.51
## cor(cmpgnMaxDisbursDvsR_mn,medianAgeYr_mn) 0.00 0.27 -0.51
## cor(primaryTotDvsR_mn,medianAgeYr_mn) 0.00 0.26 -0.50
## cor(log_popDensPerSqMi_mn,medianAgeYr_mn) 0.02 0.27 -0.50
## cor(Intercept,pctEduGradProAge25plus_mn) 0.01 0.27 -0.51
## cor(incumbentD_mn,pctEduGradProAge25plus_mn) 0.00 0.27 -0.51
## cor(incumbentR_mn,pctEduGradProAge25plus_mn) -0.00 0.27 -0.54
## cor(primaryUnopD_mn,pctEduGradProAge25plus_mn) 0.01 0.27 -0.50
## cor(primaryUnopR_mn,pctEduGradProAge25plus_mn) -0.00 0.27 -0.52
## cor(cmpgnMaxDisbursDvsR_mn,pctEduGradProAge25plus_mn) 0.00 0.27 -0.50
## cor(primaryTotDvsR_mn,pctEduGradProAge25plus_mn) -0.01 0.27 -0.51
## cor(log_popDensPerSqMi_mn,pctEduGradProAge25plus_mn) -0.04 0.27 -0.55
## cor(medianAgeYr_mn,pctEduGradProAge25plus_mn) -0.01 0.26 -0.51
## cor(Intercept,pctInLaborForceUnemp_mn) 0.00 0.27 -0.50
## cor(incumbentD_mn,pctInLaborForceUnemp_mn) -0.02 0.27 -0.54
## cor(incumbentR_mn,pctInLaborForceUnemp_mn) 0.01 0.27 -0.50
## cor(primaryUnopD_mn,pctInLaborForceUnemp_mn) -0.02 0.27 -0.53
## cor(primaryUnopR_mn,pctInLaborForceUnemp_mn) 0.00 0.26 -0.51
## cor(cmpgnMaxDisbursDvsR_mn,pctInLaborForceUnemp_mn) -0.01 0.26 -0.51
## cor(primaryTotDvsR_mn,pctInLaborForceUnemp_mn) 0.00 0.26 -0.50
## cor(log_popDensPerSqMi_mn,pctInLaborForceUnemp_mn) -0.01 0.27 -0.53
## cor(medianAgeYr_mn,pctInLaborForceUnemp_mn) 0.00 0.26 -0.50
## cor(pctEduGradProAge25plus_mn,pctInLaborForceUnemp_mn) 0.02 0.26 -0.48
## cor(Intercept,perCapitaIncome_mn) 0.01 0.27 -0.51
## cor(incumbentD_mn,perCapitaIncome_mn) 0.01 0.26 -0.49
## cor(incumbentR_mn,perCapitaIncome_mn) -0.02 0.26 -0.51
## cor(primaryUnopD_mn,perCapitaIncome_mn) 0.02 0.26 -0.49
## cor(primaryUnopR_mn,perCapitaIncome_mn) 0.01 0.27 -0.50
## cor(cmpgnMaxDisbursDvsR_mn,perCapitaIncome_mn) 0.00 0.27 -0.51
## cor(primaryTotDvsR_mn,perCapitaIncome_mn) 0.01 0.27 -0.51
## cor(log_popDensPerSqMi_mn,perCapitaIncome_mn) -0.03 0.27 -0.53
## cor(medianAgeYr_mn,perCapitaIncome_mn) 0.01 0.26 -0.50
## cor(pctEduGradProAge25plus_mn,perCapitaIncome_mn) -0.05 0.27 -0.55
## cor(pctInLaborForceUnemp_mn,perCapitaIncome_mn) 0.02 0.27 -0.49
## cor(Intercept,pctWhiteAlone_mn) -0.01 0.27 -0.52
## cor(incumbentD_mn,pctWhiteAlone_mn) 0.00 0.26 -0.50
## cor(incumbentR_mn,pctWhiteAlone_mn) -0.00 0.26 -0.51
## cor(primaryUnopD_mn,pctWhiteAlone_mn) 0.03 0.26 -0.48
## cor(primaryUnopR_mn,pctWhiteAlone_mn) -0.00 0.27 -0.52
## cor(cmpgnMaxDisbursDvsR_mn,pctWhiteAlone_mn) -0.00 0.27 -0.51
## cor(primaryTotDvsR_mn,pctWhiteAlone_mn) -0.01 0.26 -0.51
## cor(log_popDensPerSqMi_mn,pctWhiteAlone_mn) 0.07 0.27 -0.46
## cor(medianAgeYr_mn,pctWhiteAlone_mn) -0.01 0.27 -0.53
## cor(pctEduGradProAge25plus_mn,pctWhiteAlone_mn) 0.01 0.26 -0.51
## cor(pctInLaborForceUnemp_mn,pctWhiteAlone_mn) 0.01 0.26 -0.50
## cor(perCapitaIncome_mn,pctWhiteAlone_mn) 0.02 0.27 -0.50
## u-95% CI Eff.Sample Rhat
## sd(Intercept) 3.93 2457 1.00
## sd(incumbentD_mn) 7.11 3072 1.00
## sd(incumbentR_mn) 7.50 1755 1.00
## sd(primaryUnopD_mn) 10.25 2216 1.00
## sd(primaryUnopR_mn) 7.48 3515 1.00
## sd(cmpgnMaxDisbursDvsR_mn) 5.13 3781 1.00
## sd(primaryTotDvsR_mn) 6.33 2312 1.00
## sd(log_popDensPerSqMi_mn) 3.92 2625 1.00
## sd(medianAgeYr_mn) 3.93 2544 1.00
## sd(pctEduGradProAge25plus_mn) 3.25 4076 1.00
## sd(pctInLaborForceUnemp_mn) 3.98 3374 1.00
## sd(perCapitaIncome_mn) 3.64 3134 1.00
## sd(pctWhiteAlone_mn) 4.87 2483 1.00
## cor(Intercept,incumbentD_mn) 0.49 5118 1.00
## cor(Intercept,incumbentR_mn) 0.49 4019 1.00
## cor(incumbentD_mn,incumbentR_mn) 0.50 4626 1.00
## cor(Intercept,primaryUnopD_mn) 0.52 4737 1.00
## cor(incumbentD_mn,primaryUnopD_mn) 0.51 4853 1.00
## cor(incumbentR_mn,primaryUnopD_mn) 0.53 4844 1.00
## cor(Intercept,primaryUnopR_mn) 0.51 5309 1.00
## cor(incumbentD_mn,primaryUnopR_mn) 0.49 5622 1.00
## cor(incumbentR_mn,primaryUnopR_mn) 0.50 5359 1.00
## cor(primaryUnopD_mn,primaryUnopR_mn) 0.49 5230 1.00
## cor(Intercept,cmpgnMaxDisbursDvsR_mn) 0.51 5177 1.00
## cor(incumbentD_mn,cmpgnMaxDisbursDvsR_mn) 0.50 5619 1.00
## cor(incumbentR_mn,cmpgnMaxDisbursDvsR_mn) 0.52 5748 1.00
## cor(primaryUnopD_mn,cmpgnMaxDisbursDvsR_mn) 0.48 5588 1.00
## cor(primaryUnopR_mn,cmpgnMaxDisbursDvsR_mn) 0.51 5740 1.00
## cor(Intercept,primaryTotDvsR_mn) 0.50 5478 1.00
## cor(incumbentD_mn,primaryTotDvsR_mn) 0.49 5142 1.00
## cor(incumbentR_mn,primaryTotDvsR_mn) 0.51 5354 1.00
## cor(primaryUnopD_mn,primaryTotDvsR_mn) 0.49 4755 1.00
## cor(primaryUnopR_mn,primaryTotDvsR_mn) 0.52 5217 1.00
## cor(cmpgnMaxDisbursDvsR_mn,primaryTotDvsR_mn) 0.50 5191 1.00
## cor(Intercept,log_popDensPerSqMi_mn) 0.54 4422 1.00
## cor(incumbentD_mn,log_popDensPerSqMi_mn) 0.52 4989 1.00
## cor(incumbentR_mn,log_popDensPerSqMi_mn) 0.56 5255 1.00
## cor(primaryUnopD_mn,log_popDensPerSqMi_mn) 0.55 4612 1.00
## cor(primaryUnopR_mn,log_popDensPerSqMi_mn) 0.49 5039 1.00
## cor(cmpgnMaxDisbursDvsR_mn,log_popDensPerSqMi_mn) 0.49 5413 1.00
## cor(primaryTotDvsR_mn,log_popDensPerSqMi_mn) 0.48 4797 1.00
## cor(Intercept,medianAgeYr_mn) 0.52 4803 1.00
## cor(incumbentD_mn,medianAgeYr_mn) 0.52 4789 1.00
## cor(incumbentR_mn,medianAgeYr_mn) 0.49 4944 1.00
## cor(primaryUnopD_mn,medianAgeYr_mn) 0.51 5031 1.00
## cor(primaryUnopR_mn,medianAgeYr_mn) 0.50 5096 1.00
## cor(cmpgnMaxDisbursDvsR_mn,medianAgeYr_mn) 0.51 4976 1.00
## cor(primaryTotDvsR_mn,medianAgeYr_mn) 0.53 5051 1.00
## cor(log_popDensPerSqMi_mn,medianAgeYr_mn) 0.54 5068 1.00
## cor(Intercept,pctEduGradProAge25plus_mn) 0.53 5657 1.00
## cor(incumbentD_mn,pctEduGradProAge25plus_mn) 0.52 5236 1.00
## cor(incumbentR_mn,pctEduGradProAge25plus_mn) 0.53 5689 1.00
## cor(primaryUnopD_mn,pctEduGradProAge25plus_mn) 0.53 5553 1.00
## cor(primaryUnopR_mn,pctEduGradProAge25plus_mn) 0.53 5134 1.00
## cor(cmpgnMaxDisbursDvsR_mn,pctEduGradProAge25plus_mn) 0.51 5584 1.00
## cor(primaryTotDvsR_mn,pctEduGradProAge25plus_mn) 0.51 5445 1.00
## cor(log_popDensPerSqMi_mn,pctEduGradProAge25plus_mn) 0.48 5583 1.00
## cor(medianAgeYr_mn,pctEduGradProAge25plus_mn) 0.52 5507 1.00
## cor(Intercept,pctInLaborForceUnemp_mn) 0.52 5552 1.00
## cor(incumbentD_mn,pctInLaborForceUnemp_mn) 0.50 5477 1.00
## cor(incumbentR_mn,pctInLaborForceUnemp_mn) 0.52 5537 1.00
## cor(primaryUnopD_mn,pctInLaborForceUnemp_mn) 0.51 5358 1.00
## cor(primaryUnopR_mn,pctInLaborForceUnemp_mn) 0.51 5504 1.00
## cor(cmpgnMaxDisbursDvsR_mn,pctInLaborForceUnemp_mn) 0.50 4946 1.00
## cor(primaryTotDvsR_mn,pctInLaborForceUnemp_mn) 0.50 5435 1.00
## cor(log_popDensPerSqMi_mn,pctInLaborForceUnemp_mn) 0.51 4900 1.00
## cor(medianAgeYr_mn,pctInLaborForceUnemp_mn) 0.51 5177 1.00
## cor(pctEduGradProAge25plus_mn,pctInLaborForceUnemp_mn) 0.51 5447 1.00
## cor(Intercept,perCapitaIncome_mn) 0.52 5678 1.00
## cor(incumbentD_mn,perCapitaIncome_mn) 0.52 5070 1.00
## cor(incumbentR_mn,perCapitaIncome_mn) 0.51 5616 1.00
## cor(primaryUnopD_mn,perCapitaIncome_mn) 0.52 5637 1.00
## cor(primaryUnopR_mn,perCapitaIncome_mn) 0.54 5607 1.00
## cor(cmpgnMaxDisbursDvsR_mn,perCapitaIncome_mn) 0.52 5263 1.00
## cor(primaryTotDvsR_mn,perCapitaIncome_mn) 0.51 5331 1.00
## cor(log_popDensPerSqMi_mn,perCapitaIncome_mn) 0.49 5618 1.00
## cor(medianAgeYr_mn,perCapitaIncome_mn) 0.52 5336 1.00
## cor(pctEduGradProAge25plus_mn,perCapitaIncome_mn) 0.48 5108 1.00
## cor(pctInLaborForceUnemp_mn,perCapitaIncome_mn) 0.53 5449 1.00
## cor(Intercept,pctWhiteAlone_mn) 0.52 4854 1.00
## cor(incumbentD_mn,pctWhiteAlone_mn) 0.52 4475 1.00
## cor(incumbentR_mn,pctWhiteAlone_mn) 0.50 5045 1.00
## cor(primaryUnopD_mn,pctWhiteAlone_mn) 0.53 5200 1.00
## cor(primaryUnopR_mn,pctWhiteAlone_mn) 0.51 5270 1.00
## cor(cmpgnMaxDisbursDvsR_mn,pctWhiteAlone_mn) 0.51 5504 1.00
## cor(primaryTotDvsR_mn,pctWhiteAlone_mn) 0.50 5560 1.00
## cor(log_popDensPerSqMi_mn,pctWhiteAlone_mn) 0.58 4314 1.00
## cor(medianAgeYr_mn,pctWhiteAlone_mn) 0.50 5622 1.00
## cor(pctEduGradProAge25plus_mn,pctWhiteAlone_mn) 0.51 5682 1.00
## cor(pctInLaborForceUnemp_mn,pctWhiteAlone_mn) 0.53 5641 1.00
## cor(perCapitaIncome_mn,pctWhiteAlone_mn) 0.53 4954 1.00
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept -1.02 2.94 -6.89 4.93 4401 1.00
## incumbentD_mn 13.85 4.41 6.38 23.50 3511 1.00
## incumbentR_mn -7.83 4.36 -17.29 -0.23 3488 1.00
## primaryUnopD_mn 1.23 3.28 -5.25 7.86 3487 1.00
## primaryUnopR_mn -2.25 3.46 -9.48 4.32 3545 1.00
## cmpgnMaxDisbursDvsR_mn 5.06 2.13 1.61 10.02 3508 1.00
## primaryTotDvsR_mn 4.30 2.35 0.35 9.55 2798 1.00
## log_popDensPerSqMi_mn 1.01 1.16 -1.10 3.48 3773 1.00
## medianAgeYr_mn -0.17 0.93 -2.06 1.67 3662 1.00
## pctEduGradProAge25plus_mn -0.94 2.08 -5.10 3.26 3605 1.00
## pctInLaborForceUnemp_mn 1.10 1.50 -1.72 4.27 3759 1.00
## perCapitaIncome_mn 2.54 2.26 -1.74 7.26 3213 1.00
## pctWhiteAlone_mn -2.22 1.67 -5.84 0.75 3717 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).The “less regularized” varying-effects model appears to have very similar effective sample sizes, sd, and cor parameter estimates compared to the “more regularized” model, and again the R-hat values (and posterior density overlay / MCMC chain trace plots) indicate a good model fit.
Let’s compare the population-level parameter posterior estimates between the two varying-effects models:
While the fixed-effects models were quite similar between the greater and lesser regularization priors, the varying-effects model with less regularization shows a much greater polarization of posterior estimates (negative estimates and positive estimates extended further in those respective directions), and also much wider uncertainty (with 90% credible intervals that mostly entirely encompass those of the “more regularization” model).
In short, the fixed-effects model does not seem nearly as influenced by choices in prior regularization as the varying-effects model.
More comparisons among No-interactions models
The next line of code evaluates the effective number of parameters for each model. We’d expect the greatest number for the varying-effects models, since there are varying intercepts and slopes.
While I’d like to use PSIS-LOO for evaluating the effective number of parameters for each model, the varying-effects models take a considerably long time to run, so I’ll fall back on WAIC instead.
The plot below summarizes estimates of \(p_{WAIC}\), the estimated effective number of parameters, for each model.
For reference, there are 13 population parameters (Intercept + 13 predictors), and a theoretical maximum of 5,655 group-level effects parameters (13 * 435 [the number of districts]), without considering the between-groups/-effects correlation and standard deviation parameters. In other words, the varying-effects models could be massively overfit.
pWAIC_est <-
tibble(model = c("m2.allFull_fIfS - more regularized",
"m2.allFull_fIfS2 - less regularized",
"m2.allFull_vIvS - more regularized",
"m2.allFull_vIvS2 - less regularized"),
pWAIC = c(waic(m2.allFull_fIfS)$estimates["p_waic","Estimate"],
waic(m2.allFull_fIfS2)$estimates["p_waic","Estimate"],
waic(m2.allFull_vIvS)$estimates["p_waic","Estimate"],
waic(m2.allFull_vIvS2)$estimates["p_waic","Estimate"]),
pW_SE = c(waic(m2.allFull_fIfS)$estimates["p_waic","SE"],
waic(m2.allFull_fIfS2)$estimates["p_waic","SE"],
waic(m2.allFull_vIvS)$estimates["p_waic","SE"],
waic(m2.allFull_vIvS2)$estimates["p_waic","SE"]) )
pWAIC_est %>%
ggplot(aes(x = model, y = pWAIC,
ymin = pWAIC - 2 * pW_SE, ymax = pWAIC + 2 * pW_SE,
color = model)) +
geom_pointrange() +
coord_flip() +
labs(title = "Estimated effective # parameters \n for no-interaction models",
subtitle = expression("Point estimate \u00B1 2 standard errors"),
x = "Model", y = expression("Estimated effective # parameters - "*p[WAIC])) +
scale_color_manual(values = color_set8[c(2, 1, 4, 3)], guide = F) +
theme_ds1()
This plot seems truly impressive - while the less-regularized varying-effects model m2.allFull_vIvS2 plausibly has more effective parameters than either of the fixed-effects models, the (approximately) 95% credible intervals largely overlap, and the varying-effects models do not have so many more estimated effective parameters that it is immediately evident there is major overfitting.
Let’s compare relative model weights among these four models to better understand the relative estimated deviance in out-of-sample predictions.
(m2.allFull2_waic <- waic(m2.allFull_fIfS, m2.allFull_fIfS2,
m2.allFull_vIvS, m2.allFull_vIvS2,
compare = T) )## WAIC SE
## m2.allFull_fIfS 157.33 26.93
## m2.allFull_fIfS2 162.33 29.52
## m2.allFull_vIvS 139.16 22.81
## m2.allFull_vIvS2 118.75 20.01
## m2.allFull_fIfS - m2.allFull_fIfS2 -5.00 3.81
## m2.allFull_fIfS - m2.allFull_vIvS 18.17 9.25
## m2.allFull_fIfS - m2.allFull_vIvS2 38.58 14.42
## m2.allFull_fIfS2 - m2.allFull_vIvS 23.17 11.19
## m2.allFull_fIfS2 - m2.allFull_vIvS2 43.58 16.45
## m2.allFull_vIvS - m2.allFull_vIvS2 20.41 7.49While the more-regularized (original) fixed-effects model, m2.allFull_vIvS, has a lower (“better”) WAIC value than its less-regularized alternative, the opposite is true for the varying-effects model. There, the less-regularized (revised) model, m2.allFull_vIvS2 has a much lower WAIC value, indicating it is estimated to have less prediction deviance for out-of-sample data.
After discussion with Dr. Frey, seems very, very plausible the “less regularized” model may actually just be overfit - however, because the estimated effective number of parameters are not excessive in any case, I’ll move forward with the model(s) assigned the most WAIC value-based weighting in the next section.
This newer varying-effects model also outperforms the two fixed-effects models by WAIC value - I suspect this means it will be allocated a majority of WAIC value-based weighting.
(m2.allFull_modWts2 <- model_weights(m2.allFull_fIfS, m2.allFull_fIfS2,
m2.allFull_vIvS, m2.allFull_vIvS2,
weights = "waic") )## m2.allFull_fIfS m2.allFull_fIfS2 m2.allFull_vIvS m2.allFull_vIvS2
## 4.188e-09 3.443e-10 3.698e-05 1.000e+00The “less regularized” varying-effects model is essentially assigned all the model weight when compared with the two fixed-effects models and the original varying-effects model. Combined with the effective number of parameters estimates above, I’ll go forward with the less-regularized varying-effects model.
Quick comparison between varying-effects models
It seems clear that the less-strict priors of the revised m2.allFull_vIvS2 is estimated to result in improved out-of-sample estimates, but how does this relate to in-sample inference?
# prediction of each district (columns) for each post-warmup iteration (rows)
postPred_m2.allFull_vIvS <- posterior_predict(m2.allFull_vIvS)
# mean for each per-year district (column of postPred_m2.allFull_vIvS)
meanPred_m2.allFull_vIvS <- apply(postPred_m2.allFull_vIvS, 2, mean)
predVsObs.2v <-
tibble(district2 = houseElectionDat_std_agg$district2,
generalWinD = houseElectionDat_std_agg$generalWinD,
meanPredWinD = meanPred_m2.allFull_vIvS )
nBins <- 100
cols <- c("darkred", "red",
"lightgrey",
"blue", "darkblue")
colGradient <- colorRampPalette(cols)
cut.cols <- colGradient(50)
cuts <- cut(predVsObs.2v$meanPredWinD, nBins)
names(cuts) <- sapply(cuts,function(t)
cut.cols[which(as.character(t) == levels(cuts))])
predVsObs.2v %>%
ggplot(aes(x = meanPredWinD, fill = cut(meanPredWinD, nBins))) +
geom_histogram(bins = nBins) +
facet_wrap(. ~ generalWinD) +
labs(title = "Posterior prediction check for m2.allFull_vIvS",
subtitle = "Faceted by # Democratic candidates elected (out of 3 elections: 2012/14/16)",
x = "Mean count of model predictions, D candidate elected",
y = "Count") +
scale_fill_manual(values = cut.cols,
guide = F) +
theme_ds1()
# prediction of each district (columns) for each post-warmup iteration (rows)
postPred_m2.allFull_vIvS2 <- posterior_predict(m2.allFull_vIvS2)
# mean for each per-year district (column of postPred_m2.allFull_vIvS2)
meanPred_m2.allFull_vIvS2 <- apply(postPred_m2.allFull_vIvS2, 2, mean)
predVsObs.2v2 <-
tibble(district2 = houseElectionDat_std_agg$district2,
generalWinD = houseElectionDat_std_agg$generalWinD,
meanPredWinD = meanPred_m2.allFull_vIvS2 )
nBins <- 100
cols <- c("darkred", "red",
"lightgrey",
"blue", "darkblue")
colGradient <- colorRampPalette(cols)
cut.cols <- colGradient(50)
cuts <- cut(predVsObs.2v2$meanPredWinD, nBins)
names(cuts) <- sapply(cuts,function(t)
cut.cols[which(as.character(t) == levels(cuts))])
predVsObs.2v2 %>%
ggplot(aes(x = meanPredWinD, fill = cut(meanPredWinD, nBins))) +
geom_histogram(bins = nBins) +
facet_wrap(. ~ generalWinD) +
labs(title = "Posterior prediction check for m2.allFull_vIvS2",
subtitle = "Faceted by # Democratic candidates elected (out of 3 elections: 2012/14/16)",
x = "Mean count of model predictions, D candidate elected",
y = "Count") +
scale_fill_manual(values = cut.cols,
guide = F) +
theme_ds1()
# comparison of fixed- and varying-effects model predictions
ggplot() +
geom_density(data = predVsObs.2v,
aes(x = meanPredWinD), fill = color_set8[3], alpha = 0.5,
adjust = 0.1) +
geom_density(data = predVsObs.2v2,
aes(x = meanPredWinD), fill = color_set8[4], alpha = 0.5,
adjust = 0.1) +
labs(title = "Comparing posterior prediction densities, varying-effects models",
subtitle = "More (Lt green) vs Less regularized priors (Dk green)",
x = "Mean proportion of model predictions, D candidate elected",
y = "Density") +
theme_ds1()
The less-regularized priors model m2.allFullvIvS2 tightens predictions very considerably around the observed count of Democratic candidates elected across the three elections, in each of the four plot panels.
Just for quantification’s sake, let’s again compare the two varying-effects models, again using the “miss” metric of a predicted count being 0.5 or more away from the observed count:
| generalWinD | hit.v | miss.v | hit.v2 | miss.v2 |
|---|---|---|---|---|
| 0 | 225 | 0 | 225 | 0 |
| 1 | 12 | 9 | 21 | 0 |
| 2 | 4 | 1 | 5 | 0 |
| 3 | 183 | 1 | 184 | 0 |
There isn’t a single observation with a model-prediction “miss” of 0.5 or more for the less-regularized varying-effects model. There is every possibility the model is simply overfit, but I’ll proceed with the less-regularized varying-effects model on the off chance the model is well fit to the data.
Considering potential interactions
Moving forward with the “less-regularized priors, varying-effects” model (m2.allFull_vIvS2), let’s consider potentially incorporating interactions into the model. I will add them one at a time, evaluate their posterior distribution estimates, and consider adding each predictor with preference to those with posterior credible intervals not concentrated about 0.
The following predictors seem to me to be of interest:
cmpgnMaxDisbursDvsR\(\times\)primaryTotDvsR- evaluating the interaction between maximum campaign-reported spending differential and primary vote totals differentialperCapitaIncome\(\times\)pctEduGradProAge25plus- evaluating the interaction between per capita income and proportion of the age 25+ population with a graduate/professional degreeperCapitaIncome\(\times\)pctInLaborForceUnemp- evaluating the interaction between per capita income and proportion of the labor-force population that is unemployedincumbentD\(\times\)cmpgnMaxDisbursDvsR,incumbentR\(\times\)cmpgnMaxDisbursDvsR- evaluating the interaction between party incumbency and maximum reported spending differentialincumbentD\(\times\)primaryTotDvsR,incumbentR\(\times\)primaryTotDvsR- evaluating the interaction between party incumbency and total primary votes differential
m2.allFull_vIvS2_int1 <-
update(m2.allFull_vIvS2, . ~ . + cmpgnMaxDisbursDvsR_mn:primaryTotDvsR_mn,
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m2.allFull_vIvS2_int1")
m2.allFull_vIvS2_int2 <-
update(m2.allFull_vIvS2, . ~ . + perCapitaIncome_mn:pctEduGradProAge25plus_mn,
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m2.allFull_vIvS2_int2")
m2.allFull_vIvS2_int3 <-
update(m2.allFull_vIvS2, . ~ . + perCapitaIncome_mn:pctInLaborForceUnemp_mn,
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m2.allFull_vIvS2_int3")
m2.allFull_vIvS2_int4 <-
update(m2.allFull_vIvS2, . ~ . + incumbentD_mn:cmpgnMaxDisbursDvsR_mn +
incumbentR_mn:cmpgnMaxDisbursDvsR_mn,
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m2.allFull_vIvS2_int4")
m2.allFull_vIvS2_int5 <-
update(m2.allFull_vIvS2, . ~ . + incumbentD_mn:primaryTotDvsR_mn +
incumbentR_mn:primaryTotDvsR_mn,
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m2.allFull_vIvS2_int5")
While most of the interactions considered appear to very plausibly have no marginal predicted association with the election outcome, the top two interactions do have relatively narrow credible intervals, and the first seems to quite plausibly have a nonzero marginal effect.
Therefore, I will retain the first interaction model, and then compare this interaction model against m2.allFull_vIvS2.
(m2.allFull_vIvS2_modWts <- model_weights(m2.allFull_vIvS2,
m2.allFull_vIvS2_int1,
weights = "waic") )## m2.allFull_vIvS2 m2.allFull_vIvS2_int1
## 0.2722 0.7278Comparing this latest model against the less-regularized varying-effects model, the newer model featuring an interaction is “assigned” roughly 73% of the relative model weight based on WAIC values.
I take the relative weights being shared to such an extent between the two models as evidence that relying entirely on the model containing the interaction effectively places too much weight on the marginal contribution of the interaction.
To avoid this overconfidence, I will generate the “final” model using model-weighted averaging between these two models, with relative weights equal to what we see above (about a 25/75 split). This will pull posterior estimates to approximately the middle of the posterior distribution of the non-interaction model and that of the interaction-including model for each parameter.
Let’s first check the summary of the new interaction model.
summary(m2.allFull_vIvS2_int1)## Family: binomial
## Links: mu = logit
## Formula: generalWinD | trials(nObs) ~ incumbentD_mn + incumbentR_mn + primaryUnopD_mn + primaryUnopR_mn + cmpgnMaxDisbursDvsR_mn + primaryTotDvsR_mn + log_popDensPerSqMi_mn + medianAgeYr_mn + pctEduGradProAge25plus_mn + pctInLaborForceUnemp_mn + perCapitaIncome_mn + pctWhiteAlone_mn + (1 + incumbentD_mn + incumbentR_mn + primaryUnopD_mn + primaryUnopR_mn + cmpgnMaxDisbursDvsR_mn + primaryTotDvsR_mn + log_popDensPerSqMi_mn + medianAgeYr_mn + pctEduGradProAge25plus_mn + pctInLaborForceUnemp_mn + perCapitaIncome_mn + pctWhiteAlone_mn | district2) + cmpgnMaxDisbursDvsR_mn:primaryTotDvsR_mn
## Data: houseElectionDat_std_agg (Number of observations: 435)
## Samples: 3 chains, each with iter = 7000; warmup = 3500; thin = 2;
## total post-warmup samples = 5250
##
## Group-Level Effects:
## ~district2 (Number of levels: 435)
## Estimate Est.Error l-95% CI
## sd(Intercept) 1.46 1.13 0.06
## sd(incumbentD_mn) 2.70 2.05 0.11
## sd(incumbentR_mn) 3.39 2.11 0.19
## sd(primaryUnopD_mn) 4.06 2.94 0.18
## sd(primaryUnopR_mn) 2.67 2.17 0.11
## sd(cmpgnMaxDisbursDvsR_mn) 1.80 1.46 0.07
## sd(primaryTotDvsR_mn) 2.33 1.77 0.10
## sd(log_popDensPerSqMi_mn) 1.39 1.11 0.05
## sd(medianAgeYr_mn) 1.59 1.11 0.07
## sd(pctEduGradProAge25plus_mn) 1.14 0.92 0.05
## sd(pctInLaborForceUnemp_mn) 1.34 1.10 0.04
## sd(perCapitaIncome_mn) 1.31 1.04 0.05
## sd(pctWhiteAlone_mn) 1.77 1.31 0.07
## cor(Intercept,incumbentD_mn) -0.03 0.27 -0.55
## cor(Intercept,incumbentR_mn) -0.04 0.27 -0.55
## cor(incumbentD_mn,incumbentR_mn) 0.00 0.27 -0.50
## cor(Intercept,primaryUnopD_mn) 0.00 0.27 -0.51
## cor(incumbentD_mn,primaryUnopD_mn) -0.01 0.27 -0.51
## cor(incumbentR_mn,primaryUnopD_mn) 0.03 0.26 -0.48
## cor(Intercept,primaryUnopR_mn) -0.02 0.27 -0.53
## cor(incumbentD_mn,primaryUnopR_mn) -0.01 0.27 -0.51
## cor(incumbentR_mn,primaryUnopR_mn) -0.03 0.27 -0.54
## cor(primaryUnopD_mn,primaryUnopR_mn) -0.05 0.27 -0.56
## cor(Intercept,cmpgnMaxDisbursDvsR_mn) -0.01 0.27 -0.52
## cor(incumbentD_mn,cmpgnMaxDisbursDvsR_mn) -0.02 0.27 -0.54
## cor(incumbentR_mn,cmpgnMaxDisbursDvsR_mn) 0.01 0.27 -0.49
## cor(primaryUnopD_mn,cmpgnMaxDisbursDvsR_mn) -0.03 0.27 -0.54
## cor(primaryUnopR_mn,cmpgnMaxDisbursDvsR_mn) -0.01 0.27 -0.52
## cor(Intercept,primaryTotDvsR_mn) -0.01 0.27 -0.53
## cor(incumbentD_mn,primaryTotDvsR_mn) -0.02 0.26 -0.51
## cor(incumbentR_mn,primaryTotDvsR_mn) -0.00 0.27 -0.52
## cor(primaryUnopD_mn,primaryTotDvsR_mn) -0.03 0.27 -0.55
## cor(primaryUnopR_mn,primaryTotDvsR_mn) -0.00 0.27 -0.52
## cor(cmpgnMaxDisbursDvsR_mn,primaryTotDvsR_mn) -0.02 0.27 -0.54
## cor(Intercept,log_popDensPerSqMi_mn) 0.04 0.27 -0.47
## cor(incumbentD_mn,log_popDensPerSqMi_mn) 0.00 0.27 -0.50
## cor(incumbentR_mn,log_popDensPerSqMi_mn) 0.06 0.27 -0.46
## cor(primaryUnopD_mn,log_popDensPerSqMi_mn) 0.04 0.27 -0.49
## cor(primaryUnopR_mn,log_popDensPerSqMi_mn) -0.01 0.27 -0.53
## cor(cmpgnMaxDisbursDvsR_mn,log_popDensPerSqMi_mn) -0.02 0.26 -0.52
## cor(primaryTotDvsR_mn,log_popDensPerSqMi_mn) -0.04 0.27 -0.55
## cor(Intercept,medianAgeYr_mn) 0.00 0.26 -0.50
## cor(incumbentD_mn,medianAgeYr_mn) 0.01 0.27 -0.50
## cor(incumbentR_mn,medianAgeYr_mn) -0.00 0.26 -0.53
## cor(primaryUnopD_mn,medianAgeYr_mn) -0.01 0.26 -0.52
## cor(primaryUnopR_mn,medianAgeYr_mn) -0.00 0.27 -0.52
## cor(cmpgnMaxDisbursDvsR_mn,medianAgeYr_mn) 0.01 0.26 -0.49
## cor(primaryTotDvsR_mn,medianAgeYr_mn) 0.01 0.26 -0.50
## cor(log_popDensPerSqMi_mn,medianAgeYr_mn) 0.02 0.26 -0.50
## cor(Intercept,pctEduGradProAge25plus_mn) 0.01 0.27 -0.52
## cor(incumbentD_mn,pctEduGradProAge25plus_mn) -0.00 0.27 -0.51
## cor(incumbentR_mn,pctEduGradProAge25plus_mn) -0.00 0.26 -0.52
## cor(primaryUnopD_mn,pctEduGradProAge25plus_mn) 0.01 0.27 -0.51
## cor(primaryUnopR_mn,pctEduGradProAge25plus_mn) -0.00 0.27 -0.51
## cor(cmpgnMaxDisbursDvsR_mn,pctEduGradProAge25plus_mn) 0.00 0.27 -0.52
## cor(primaryTotDvsR_mn,pctEduGradProAge25plus_mn) -0.01 0.27 -0.53
## cor(log_popDensPerSqMi_mn,pctEduGradProAge25plus_mn) -0.03 0.27 -0.53
## cor(medianAgeYr_mn,pctEduGradProAge25plus_mn) -0.00 0.27 -0.51
## cor(Intercept,pctInLaborForceUnemp_mn) -0.01 0.26 -0.53
## cor(incumbentD_mn,pctInLaborForceUnemp_mn) -0.02 0.27 -0.52
## cor(incumbentR_mn,pctInLaborForceUnemp_mn) 0.02 0.26 -0.49
## cor(primaryUnopD_mn,pctInLaborForceUnemp_mn) -0.01 0.27 -0.52
## cor(primaryUnopR_mn,pctInLaborForceUnemp_mn) 0.01 0.27 -0.51
## cor(cmpgnMaxDisbursDvsR_mn,pctInLaborForceUnemp_mn) -0.01 0.27 -0.52
## cor(primaryTotDvsR_mn,pctInLaborForceUnemp_mn) -0.00 0.26 -0.51
## cor(log_popDensPerSqMi_mn,pctInLaborForceUnemp_mn) -0.01 0.26 -0.52
## cor(medianAgeYr_mn,pctInLaborForceUnemp_mn) -0.00 0.27 -0.52
## cor(pctEduGradProAge25plus_mn,pctInLaborForceUnemp_mn) 0.01 0.27 -0.50
## cor(Intercept,perCapitaIncome_mn) 0.01 0.26 -0.49
## cor(incumbentD_mn,perCapitaIncome_mn) 0.01 0.27 -0.50
## cor(incumbentR_mn,perCapitaIncome_mn) -0.02 0.27 -0.53
## cor(primaryUnopD_mn,perCapitaIncome_mn) 0.01 0.27 -0.50
## cor(primaryUnopR_mn,perCapitaIncome_mn) 0.01 0.27 -0.52
## cor(cmpgnMaxDisbursDvsR_mn,perCapitaIncome_mn) 0.00 0.26 -0.49
## cor(primaryTotDvsR_mn,perCapitaIncome_mn) -0.00 0.27 -0.52
## cor(log_popDensPerSqMi_mn,perCapitaIncome_mn) -0.03 0.27 -0.54
## cor(medianAgeYr_mn,perCapitaIncome_mn) 0.00 0.26 -0.51
## cor(pctEduGradProAge25plus_mn,perCapitaIncome_mn) -0.04 0.27 -0.54
## cor(pctInLaborForceUnemp_mn,perCapitaIncome_mn) 0.02 0.27 -0.50
## cor(Intercept,pctWhiteAlone_mn) 0.00 0.26 -0.50
## cor(incumbentD_mn,pctWhiteAlone_mn) 0.01 0.27 -0.50
## cor(incumbentR_mn,pctWhiteAlone_mn) 0.00 0.27 -0.51
## cor(primaryUnopD_mn,pctWhiteAlone_mn) 0.02 0.27 -0.49
## cor(primaryUnopR_mn,pctWhiteAlone_mn) -0.00 0.27 -0.51
## cor(cmpgnMaxDisbursDvsR_mn,pctWhiteAlone_mn) -0.00 0.27 -0.52
## cor(primaryTotDvsR_mn,pctWhiteAlone_mn) 0.00 0.26 -0.51
## cor(log_popDensPerSqMi_mn,pctWhiteAlone_mn) 0.07 0.27 -0.48
## cor(medianAgeYr_mn,pctWhiteAlone_mn) -0.00 0.27 -0.53
## cor(pctEduGradProAge25plus_mn,pctWhiteAlone_mn) 0.02 0.27 -0.50
## cor(pctInLaborForceUnemp_mn,pctWhiteAlone_mn) 0.01 0.27 -0.50
## cor(perCapitaIncome_mn,pctWhiteAlone_mn) 0.02 0.27 -0.51
## u-95% CI Eff.Sample Rhat
## sd(Intercept) 4.23 3083 1.00
## sd(incumbentD_mn) 7.59 3707 1.00
## sd(incumbentR_mn) 8.18 2170 1.00
## sd(primaryUnopD_mn) 11.05 2609 1.00
## sd(primaryUnopR_mn) 7.88 3934 1.00
## sd(cmpgnMaxDisbursDvsR_mn) 5.62 3930 1.00
## sd(primaryTotDvsR_mn) 6.46 2979 1.00
## sd(log_popDensPerSqMi_mn) 4.21 3391 1.00
## sd(medianAgeYr_mn) 4.14 2392 1.00
## sd(pctEduGradProAge25plus_mn) 3.32 4173 1.00
## sd(pctInLaborForceUnemp_mn) 4.08 3591 1.00
## sd(perCapitaIncome_mn) 3.87 3486 1.00
## sd(pctWhiteAlone_mn) 5.00 3181 1.00
## cor(Intercept,incumbentD_mn) 0.50 4295 1.00
## cor(Intercept,incumbentR_mn) 0.49 4273 1.00
## cor(incumbentD_mn,incumbentR_mn) 0.51 4361 1.00
## cor(Intercept,primaryUnopD_mn) 0.51 4512 1.00
## cor(incumbentD_mn,primaryUnopD_mn) 0.50 4556 1.00
## cor(incumbentR_mn,primaryUnopD_mn) 0.54 4120 1.00
## cor(Intercept,primaryUnopR_mn) 0.51 4265 1.00
## cor(incumbentD_mn,primaryUnopR_mn) 0.51 3945 1.00
## cor(incumbentR_mn,primaryUnopR_mn) 0.49 4251 1.00
## cor(primaryUnopD_mn,primaryUnopR_mn) 0.49 4122 1.00
## cor(Intercept,cmpgnMaxDisbursDvsR_mn) 0.50 4519 1.00
## cor(incumbentD_mn,cmpgnMaxDisbursDvsR_mn) 0.48 4521 1.00
## cor(incumbentR_mn,cmpgnMaxDisbursDvsR_mn) 0.52 4510 1.00
## cor(primaryUnopD_mn,cmpgnMaxDisbursDvsR_mn) 0.50 4417 1.00
## cor(primaryUnopR_mn,cmpgnMaxDisbursDvsR_mn) 0.50 4163 1.00
## cor(Intercept,primaryTotDvsR_mn) 0.50 4635 1.00
## cor(incumbentD_mn,primaryTotDvsR_mn) 0.50 4655 1.00
## cor(incumbentR_mn,primaryTotDvsR_mn) 0.51 4618 1.00
## cor(primaryUnopD_mn,primaryTotDvsR_mn) 0.49 4596 1.00
## cor(primaryUnopR_mn,primaryTotDvsR_mn) 0.50 4649 1.00
## cor(cmpgnMaxDisbursDvsR_mn,primaryTotDvsR_mn) 0.50 4402 1.00
## cor(Intercept,log_popDensPerSqMi_mn) 0.54 4263 1.00
## cor(incumbentD_mn,log_popDensPerSqMi_mn) 0.52 4452 1.00
## cor(incumbentR_mn,log_popDensPerSqMi_mn) 0.55 4582 1.00
## cor(primaryUnopD_mn,log_popDensPerSqMi_mn) 0.54 4414 1.00
## cor(primaryUnopR_mn,log_popDensPerSqMi_mn) 0.49 4559 1.00
## cor(cmpgnMaxDisbursDvsR_mn,log_popDensPerSqMi_mn) 0.50 4489 1.00
## cor(primaryTotDvsR_mn,log_popDensPerSqMi_mn) 0.49 4042 1.00
## cor(Intercept,medianAgeYr_mn) 0.51 4503 1.00
## cor(incumbentD_mn,medianAgeYr_mn) 0.52 4673 1.00
## cor(incumbentR_mn,medianAgeYr_mn) 0.49 4588 1.00
## cor(primaryUnopD_mn,medianAgeYr_mn) 0.50 4823 1.00
## cor(primaryUnopR_mn,medianAgeYr_mn) 0.51 4595 1.00
## cor(cmpgnMaxDisbursDvsR_mn,medianAgeYr_mn) 0.53 4319 1.00
## cor(primaryTotDvsR_mn,medianAgeYr_mn) 0.51 4163 1.00
## cor(log_popDensPerSqMi_mn,medianAgeYr_mn) 0.52 4739 1.00
## cor(Intercept,pctEduGradProAge25plus_mn) 0.52 4034 1.00
## cor(incumbentD_mn,pctEduGradProAge25plus_mn) 0.51 4603 1.00
## cor(incumbentR_mn,pctEduGradProAge25plus_mn) 0.50 4568 1.00
## cor(primaryUnopD_mn,pctEduGradProAge25plus_mn) 0.53 4303 1.00
## cor(primaryUnopR_mn,pctEduGradProAge25plus_mn) 0.51 4437 1.00
## cor(cmpgnMaxDisbursDvsR_mn,pctEduGradProAge25plus_mn) 0.52 4192 1.00
## cor(primaryTotDvsR_mn,pctEduGradProAge25plus_mn) 0.50 4680 1.00
## cor(log_popDensPerSqMi_mn,pctEduGradProAge25plus_mn) 0.50 4444 1.00
## cor(medianAgeYr_mn,pctEduGradProAge25plus_mn) 0.51 4246 1.00
## cor(Intercept,pctInLaborForceUnemp_mn) 0.51 4143 1.00
## cor(incumbentD_mn,pctInLaborForceUnemp_mn) 0.50 4076 1.00
## cor(incumbentR_mn,pctInLaborForceUnemp_mn) 0.52 4724 1.00
## cor(primaryUnopD_mn,pctInLaborForceUnemp_mn) 0.50 4258 1.00
## cor(primaryUnopR_mn,pctInLaborForceUnemp_mn) 0.52 4374 1.00
## cor(cmpgnMaxDisbursDvsR_mn,pctInLaborForceUnemp_mn) 0.52 4578 1.00
## cor(primaryTotDvsR_mn,pctInLaborForceUnemp_mn) 0.51 4408 1.00
## cor(log_popDensPerSqMi_mn,pctInLaborForceUnemp_mn) 0.50 3744 1.00
## cor(medianAgeYr_mn,pctInLaborForceUnemp_mn) 0.51 4361 1.00
## cor(pctEduGradProAge25plus_mn,pctInLaborForceUnemp_mn) 0.53 3962 1.00
## cor(Intercept,perCapitaIncome_mn) 0.51 4082 1.00
## cor(incumbentD_mn,perCapitaIncome_mn) 0.52 4154 1.00
## cor(incumbentR_mn,perCapitaIncome_mn) 0.50 4526 1.00
## cor(primaryUnopD_mn,perCapitaIncome_mn) 0.53 4094 1.00
## cor(primaryUnopR_mn,perCapitaIncome_mn) 0.51 4232 1.00
## cor(cmpgnMaxDisbursDvsR_mn,perCapitaIncome_mn) 0.51 4529 1.00
## cor(primaryTotDvsR_mn,perCapitaIncome_mn) 0.51 4414 1.00
## cor(log_popDensPerSqMi_mn,perCapitaIncome_mn) 0.48 4649 1.00
## cor(medianAgeYr_mn,perCapitaIncome_mn) 0.51 4521 1.00
## cor(pctEduGradProAge25plus_mn,perCapitaIncome_mn) 0.49 4496 1.00
## cor(pctInLaborForceUnemp_mn,perCapitaIncome_mn) 0.53 4099 1.00
## cor(Intercept,pctWhiteAlone_mn) 0.51 4233 1.00
## cor(incumbentD_mn,pctWhiteAlone_mn) 0.53 4549 1.00
## cor(incumbentR_mn,pctWhiteAlone_mn) 0.52 4662 1.00
## cor(primaryUnopD_mn,pctWhiteAlone_mn) 0.53 4281 1.00
## cor(primaryUnopR_mn,pctWhiteAlone_mn) 0.51 4539 1.00
## cor(cmpgnMaxDisbursDvsR_mn,pctWhiteAlone_mn) 0.51 4652 1.00
## cor(primaryTotDvsR_mn,pctWhiteAlone_mn) 0.51 4567 1.00
## cor(log_popDensPerSqMi_mn,pctWhiteAlone_mn) 0.58 4495 1.00
## cor(medianAgeYr_mn,pctWhiteAlone_mn) 0.50 4361 1.00
## cor(pctEduGradProAge25plus_mn,pctWhiteAlone_mn) 0.53 4568 1.00
## cor(pctInLaborForceUnemp_mn,pctWhiteAlone_mn) 0.51 4774 1.00
## cor(perCapitaIncome_mn,pctWhiteAlone_mn) 0.52 4122 1.00
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept -0.98 3.07 -7.16 5.02 4055
## incumbentD_mn 14.31 4.52 6.60 24.14 3737
## incumbentR_mn -8.15 4.49 -17.79 -0.31 3832
## primaryUnopD_mn 1.34 3.37 -5.26 8.03 4183
## primaryUnopR_mn -2.58 3.64 -10.03 4.43 3963
## cmpgnMaxDisbursDvsR_mn 5.56 2.28 1.95 10.77 3527
## primaryTotDvsR_mn 4.69 2.46 0.62 10.35 3254
## log_popDensPerSqMi_mn 1.02 1.22 -1.23 3.69 4212
## medianAgeYr_mn -0.16 1.00 -2.24 1.78 3786
## pctEduGradProAge25plus_mn -0.80 2.15 -4.98 3.60 3929
## pctInLaborForceUnemp_mn 1.12 1.59 -1.89 4.40 4469
## perCapitaIncome_mn 2.52 2.28 -1.97 7.12 3651
## pctWhiteAlone_mn -2.41 1.81 -6.29 0.80 4187
## cmpgnMaxDisbursDvsR_mn:primaryTotDvsR_mn -0.03 2.73 -5.71 4.89 4279
## Rhat
## Intercept 1.00
## incumbentD_mn 1.00
## incumbentR_mn 1.00
## primaryUnopD_mn 1.00
## primaryUnopR_mn 1.00
## cmpgnMaxDisbursDvsR_mn 1.00
## primaryTotDvsR_mn 1.00
## log_popDensPerSqMi_mn 1.00
## medianAgeYr_mn 1.00
## pctEduGradProAge25plus_mn 1.00
## pctInLaborForceUnemp_mn 1.00
## perCapitaIncome_mn 1.00
## pctWhiteAlone_mn 1.00
## cmpgnMaxDisbursDvsR_mn:primaryTotDvsR_mn 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).While model fits (based on effective sample sizes and R-hat values) look good, parameter estimates don’t seem to have changed much with the interaction model compared to the non-interaction m2.allFull_vIvS2, especially for the sd and b_ parameters.
For example, here are the “population effect” b_ parameters for each model. The interaction effect is fixed at zero for the non-interaction model:
If the interaction effect really does have a non-zero association with the outcome generalWinD, it appears to be very dependent on the particular district(s) involved, as the population-level effect is very plausibly zero. The other population-level effects do not appear greatly changed.
Let’s check the estimated number effective number of parameters for each model, to ensure adding the interaction effect doesn’t induce any worrying increase in overfitting:
pWAIC_est <-
tibble(model = c("m2.allFull_vIvS2",
"m2.allFull_vIvS2_int1"),
pWAIC = c(waic(m2.allFull_vIvS2)$estimates["p_waic","Estimate"],
waic(m2.allFull_vIvS2_int1)$estimates["p_waic","Estimate"]),
pW_SE = c(waic(m2.allFull_vIvS2)$estimates["p_waic","SE"],
waic(m2.allFull_vIvS2_int1)$estimates["p_waic","SE"]) )
pWAIC_est %>%
ggplot(aes(x = model, y = pWAIC,
ymin = pWAIC - 2 * pW_SE, ymax = pWAIC + 2 * pW_SE,
color = model)) +
geom_pointrange() +
coord_flip() +
labs(title = "Estimated effective # parameters",
subtitle = expression("Point estimate \u00B1 2 standard errors"),
x = "Model", y = expression("Estimated effective # parameters - "*p[WAIC])) +
scale_color_manual(values = color_set8[c(3, 6)], guide = F) +
theme_ds1()
Looks good! There isn’t an appreciable increase in the estimated effective number of parameters after adding the interaction.
Now, let’s see how model-weighted averaging influences (population) posterior estimates (applicable only to parameters shared between the two models).
I believe model-weighted averaging is removing the interaction effect due to the non-interaction model not including that term, rather than proportionally down-weighting it towards zero (which I would have hoped for) - the original version of the plot below (and related data checking) indicates it isn’t present. In order to include this parameter, I refit the “no interaction” model with the interaction term having an extremely regularizing prior (Normal(0, 0.00001)) in the hopes of including the term in model-weighted averaging while still minimalizing the term in the “no-interaction” model - probably not the wisest course of action here, but c’est la vie.
# fitting a version of the "no interaction" model with an extremely skeptical
# prior on the interaction effect, to include it in mod-wtd-averaging
m2.allFull_vIvS2_v2 <-
update(m2.allFull_vIvS2, . ~ . + cmpgnMaxDisbursDvsR_mn:primaryTotDvsR_mn,
# set very strict prior to effectively assign 0 to the interaction
prior = c(set_prior("normal(0, 0.00001)", class = "b",
coef = "cmpgnMaxDisbursDvsR_mn:primaryTotDvsR_mn") ),
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m2.allFull_vIvS2_v2")
# note: WAIC-based model weights comes down to roughly
# 25% "no-interaction model" / 75% "interaction model", which is very very close
# to the original 23% "no-interaction" / 77% "interaction"
# so just sticking with the previously-run weighting
m2.final_vIvS2 <-
posterior_average(m2.allFull_vIvS2_v2, m2.allFull_vIvS2_int1,
weights = m2.allFull_vIvS2_modWts)
post_b_modWtAvg <-
posterior_summary(m2.final_vIvS2,
probs = c(0.05, 0.95))[grepl("b_",
rownames(posterior_summary(m2.final_vIvS2))),] %>%
as_tibble() %>%
mutate(par = rownames(posterior_summary(m2.final_vIvS2)[grepl("b_",
rownames(posterior_summary(m2.final_vIvS2))),]),
model = "model-weighted average" ) %>%
mutate(par = str_remove_all(par, "b_"),
par = str_remove_all(par, "_mn") )
post_b_2.v2 %>%
bind_rows(post_b_modWtAvg) %>%
ggplot() +
geom_hline(yintercept = 0, color = "gray") +
geom_pointrange(aes(x = par, y = Estimate, ymin = Q5, ymax = Q95,
color = model, group = model),
position = position_dodge(width = 0.8)) +
labs(title = "Comparison of population parameters, \n individual models and model-weighted averages",
x = "Population-level parameter",
y = "Posterior estimate", color = NULL,
caption = "Mean point estimate and 90% credible intervals") +
coord_flip() +
scale_color_manual(values = color_set8[c(3, 6, 7)]) +
theme_ds1() + theme(legend.position = "top")
As expected, since the “with interaction” model doesn’t drastically change the population-level effects, the model-weighted average estimates have point estimates between that of the two models, and credible intervals are quite comparable to that of the other two models.
Now, let’s consider counterfactual analysis (model-implied prediction for theoretical values of one predictor, holding other values constant at some level), and then look at some district-level predictors for the model-weight average “model”.
Counterfactual analysis 1 - incumbency
In order to better understand how the model behaves for changes to a given predictor, it seems reasonable to choose a few predictors over which we’ll observe model-predicted counts of democratic candidates elected (out of three elections), while holding all predictors constant except for one.
I’m not certain , but it appears that Bayesian Model Averaging is not directly possible in this case - I’m using only the interaction-including model here.
This first counterfactual plot imagines each district somehow had an incumbent from the Democratic party (left column), the Republican party (right column), or neither party (middle column) for each of the three elections, and then predicts how many elections of the three would elect a Democratic candidate - the alternative is electing a Republican candidate (since no observed data involved an Other-party election). The mean and 2 standard deviations of each district’s predictions are plotted, with lower and upper bounds at 0 and 3, respectively.
This plot is further organized into how many Democratic candidates actually were elected as the rows, and observations are color-coded by the true number of incumbents from each party across the three elections.
Note: Because post-2010 redistricting moved some incumbents into new districts, there were some cases of multi-incumbent districts.
# create new data frames for each counterfactual plot
# note that incumbency is set to D / None / R for plot
# all other predictors (save one) set to the mean value
houseElectionDat_std_agg <-
houseElectionDat_std_agg %>%
mutate(cmpgnMaxDisbursXprimaryTotDvsR =
cmpgnMaxDisbursDvsR_mn * primaryTotDvsR_mn)
# imagining all districts somehow had a Democratic party incumbent each time
newDat_incumbentD <-
tibble(district2 = houseElectionDat_std_agg$district2,
nObs = rep(3, 435),
incumbentD_mn = rep(1, 435),
incumbentR_mn = rep(0, 435),
primaryUnopD_mn =
rep(mean(houseElectionDat_std_agg$primaryUnopD_mn), 435),
primaryUnopR_mn =
rep(mean(houseElectionDat_std_agg$primaryUnopR_mn), 435),
cmpgnMaxDisbursDvsR_mn =
rep(mean(houseElectionDat_std_agg$cmpgnMaxDisbursDvsR_mn), 435),
primaryTotDvsR_mn =
rep(mean(houseElectionDat_std_agg$primaryTotDvsR_mn), 435),
log_popDensPerSqMi_mn =
rep(mean(houseElectionDat_std_agg$log_popDensPerSqMi_mn), 435),
medianAgeYr_mn =
rep(mean(houseElectionDat_std_agg$medianAgeYr_mn), 435),
pctEduGradProAge25plus_mn =
rep(mean(houseElectionDat_std_agg$pctEduGradProAge25plus_mn), 435),
pctInLaborForceUnemp_mn =
rep(mean(houseElectionDat_std_agg$pctInLaborForceUnemp_mn), 435),
perCapitaIncome_mn =
rep(mean(houseElectionDat_std_agg$perCapitaIncome_mn), 435),
pctWhiteAlone_mn =
rep(mean(houseElectionDat_std_agg$pctWhiteAlone_mn), 435) )
# imagining all districts somehow had a Republican party incumbent each time
newDat_incumbentR <-
tibble(district2 = houseElectionDat_std_agg$district2,
nObs = rep(3, 435),
incumbentD_mn = rep(0, 435),
incumbentR_mn = rep(1, 435),
primaryUnopD_mn =
rep(mean(houseElectionDat_std_agg$primaryUnopD_mn), 435),
primaryUnopR_mn =
rep(mean(houseElectionDat_std_agg$primaryUnopR_mn), 435),
cmpgnMaxDisbursDvsR_mn =
rep(mean(houseElectionDat_std_agg$cmpgnMaxDisbursDvsR_mn), 435),
primaryTotDvsR_mn =
rep(mean(houseElectionDat_std_agg$primaryTotDvsR_mn), 435),
log_popDensPerSqMi_mn =
rep(mean(houseElectionDat_std_agg$log_popDensPerSqMi_mn), 435),
medianAgeYr_mn =
rep(mean(houseElectionDat_std_agg$medianAgeYr_mn), 435),
pctEduGradProAge25plus_mn =
rep(mean(houseElectionDat_std_agg$pctEduGradProAge25plus_mn), 435),
pctInLaborForceUnemp_mn =
rep(mean(houseElectionDat_std_agg$pctInLaborForceUnemp_mn), 435),
perCapitaIncome_mn =
rep(mean(houseElectionDat_std_agg$perCapitaIncome_mn), 435),
pctWhiteAlone_mn =
rep(mean(houseElectionDat_std_agg$pctWhiteAlone_mn), 435) )
# imagining all districts somehow had neither a Democratic party incumbent
# nor a Republican party incumbent each time
newDat_incumbentN <-
tibble(district2 = houseElectionDat_std_agg$district2,
nObs = rep(3, 435),
incumbentD_mn = rep(0, 435),
incumbentR_mn = rep(0, 435),
primaryUnopD_mn =
rep(mean(houseElectionDat_std_agg$primaryUnopD_mn), 435),
primaryUnopR_mn =
rep(mean(houseElectionDat_std_agg$primaryUnopR_mn), 435),
cmpgnMaxDisbursDvsR_mn =
rep(mean(houseElectionDat_std_agg$cmpgnMaxDisbursDvsR_mn), 435),
primaryTotDvsR_mn =
rep(mean(houseElectionDat_std_agg$primaryTotDvsR_mn), 435),
log_popDensPerSqMi_mn =
rep(mean(houseElectionDat_std_agg$log_popDensPerSqMi_mn), 435),
medianAgeYr_mn =
rep(mean(houseElectionDat_std_agg$medianAgeYr_mn), 435),
pctEduGradProAge25plus_mn =
rep(mean(houseElectionDat_std_agg$pctEduGradProAge25plus_mn), 435),
pctInLaborForceUnemp_mn =
rep(mean(houseElectionDat_std_agg$pctInLaborForceUnemp_mn), 435),
perCapitaIncome_mn =
rep(mean(houseElectionDat_std_agg$perCapitaIncome_mn), 435),
pctWhiteAlone_mn =
rep(mean(houseElectionDat_std_agg$pctWhiteAlone_mn), 435) )
predCount_incD <- posterior_predict(m2.allFull_vIvS2_int1,
newdata = newDat_incumbentD)
predCount_incR <- posterior_predict(m2.allFull_vIvS2_int1,
newdata = newDat_incumbentR)
predCount_incN <- posterior_predict(m2.allFull_vIvS2_int1,
newdata = newDat_incumbentN)
# need to take the mean of each column for estimate, then SD
predCount_Dat <-
tibble(district2 = houseElectionDat_std_agg$district2,
generalWinD = houseElectionDat_std_agg$generalWinD,
incumbentD = houseElectionDat_std_agg$incumbentD_mn,
incumbentR = houseElectionDat_std_agg$incumbentR_mn,
predCt_incD = apply(predCount_incD, 2, mean),
predSD_incD = apply(predCount_incD, 2, sd),
predCt_incR = apply(predCount_incR, 2, mean),
predSD_incR = apply(predCount_incR, 2, sd),
predCt_incN = apply(predCount_incN, 2, mean),
predSD_incN = apply(predCount_incN, 2, sd) )
# see SO user hadley's (accepted) answer
# https://stackoverflow.com/questions/25925556/gather-multiple-sets-of-columns
predCount_Dat <-
predCount_Dat %>%
gather(key, value, -district2, -generalWinD, -incumbentR, -incumbentD) %>%
extract(key, into = c("metric", "incumbentParty"),
# ^ = start of string,
# ([:alpha:]*) = first group; match upper/lower alphabetic characters
# match to end of group
# \\_ = literal underscore; separates first and second group here
# ([:alpha:]*) = second group; match upper/lower alphabetic characters
# match to end of group
regex = "^([:alpha:]*)\\_([:alpha:]*)") %>%
spread(metric, value) %>%
mutate(incumbentParty = str_remove(incumbentParty, "inc"),
# maximum of predCt - 2 SDs and 0
predMin = if_else(predCt - 2*predSD < 0, 0, predCt - 2*predSD),
# minimum of predCt + 2 SDs and 3
predMax = if_else(predCt + 2*predSD > 3, 3, predCt + 2*predSD),
incFactor = factor(
if_else(incumbentR*3 == 4, "4R",
if_else(incumbentR*3 == 3 & incumbentD*3 == 0, "3R",
if_else(incumbentR*3 == 3 & incumbentD*3 == 1, "3R 1D",
if_else(incumbentR*3 == 2 & incumbentD*3 == 0, "2R",
if_else(incumbentR*3 == 2 & incumbentD*3 == 1, "2R 1D",
if_else(incumbentR*3 == 1 & incumbentD*3 == 0, "1R",
if_else(incumbentR*3 == 0 & incumbentD*3 == 0, "None",
if_else(incumbentR*3 == 1 & incumbentD*3 == 1, "1D 1R",
if_else(incumbentR*3 == 2 & incumbentD*3 == 2, "2R 2D",
if_else(incumbentR*3 == 0 & incumbentD*3 == 1, "1D",
if_else(incumbentR*3 == 1 & incumbentD*3 == 2, "2D 1R",
if_else(incumbentR*3 == 0 & incumbentD*3 == 2, "2D",
if_else(incumbentR*3 == 1 & incumbentD*3 == 3, "3D 1R",
if_else(incumbentR*3 == 0 & incumbentD*3 == 3, "3D",
if_else(incumbentD*3 == 4, "4D", "CHECK"))))))))))))))),
levels = c("4R", "3R", "3R 1D", "2R", "2R 1D", "1R",
"None", "1D 1R", "2R 2D", "1D", "2D 1R",
"2D", "3D 1R", "3D", "4D", "CHECK"),
ordered = T) )
library(RColorBrewer)
# take the first -6- reds; last two are too light to be distinguishing
plotReds <- rev(brewer.pal(8, name = "Reds")) # original order is light to dark
# take the last two values, and use a grey for "none"
plotPurps <- brewer.pal(4, name = "Purples")
# opposite of plotReds, take the last -6- blues
plotBlues <- brewer.pal(8, name = "Blues")
# note: plot will use all 6 Reds, no Grey (for "None"),
# 2 Purples (for equal # D/R), and all 5 of 6 Blues (no 4th / #6)
plotCols <- c(plotReds[1:6], plotPurps[3:4], plotBlues[c(3:5,7,8)])
predCount_Dat %>%
ggplot(aes(x = district2, y = predCt, ymin = predMin, ymax = predMax,
color = incFactor)) +
geom_pointrange() +
facet_wrap(generalWinD ~ incumbentParty, ncol = 3) +
scale_color_manual(values = plotCols) +
#scale_color_viridis_d() +
labs(title = "Counterfactual plot considering 'all-[party] incumbents'; Pt Est \u00B1 2 SDs ",
subtitle = "Numbers: Actual # of Democratic candidates elected across the three elections \n Letters: Dem / No / Rep incumbent assigned for all 3 elections \n Colors: Actual # incumbents by party across the three elections",
x = NULL, y = "Predicted # Democratic candidates elected", color = NULL,
caption = "Due to redistricting in 2012, some districts had >3 incumbents") +
theme_ds1() +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
This is a very potentially-confusing plot, and I’ll do my best to discuss what I think it is showing.
Across each row the same districts are plotted with different counterfactual data assigned - the sample-wide mean of all predictors except for incumbency, where each district was manually assigned to have either a Democratic incumbent (left column), no incumbent (middle), or a Republican incumbent (right) heading into each election.
Within each column, districts are divided based on how many times in the three elections a Democratic candidate was actually elected (0 - 3).
This plot seems to show that the model considers incumbency for one party or the other is a very strongly associated predictor with the actual count of that party’s candidates being elected, which seems sensible. There appears to be greater relative weight toward Democratic candidate incumbency, as the “\(\pm 2 ~\text{SD}\)” prediction ranges are wider for the “assigned all-R incumbents” cells per row than for their “assigned all-D incumbents” counterparts. There is considerable uncertainty when no incumbents are assigned across the elections, and most districts seem to have a point estimate of 1.8 or so Democratic candidates being elected.
Also, there is evidence that districts are inferred to have different underlying proclivities for electing a given party’s candidates - in the middle column it is apparent that point estimates differ somewhat, while in the left and right columns one can see the ranges differ.
Counterfactual analysis 2 - difference in maximum campaign-reported spend
Note: In the “Summary table of pre-standardization values shown earlier, each year’s election’s mean difference in Democratic-minus-Republican maximum spend has Republicans with a roughly $210-240,000 spending advantage, with a standard deviation of $1.7-1.9 million.
This counterfactual plot also assigns non-incumbency predictors to sample-wide mean values, and considers predicted # of Democratic candidates elected (again, the alternative is the Republican candidate being elected) with the manipulated variable being the number of standard deviations from the population mean for difference in campaign-reported spend for the Democratic party candidate vs. the Republican party candidate. Incumbency counts are allowed to remain at the actual district-level values. And once more, rows indicate the observed number of Democratic candidates elected across the three elections.
# create new data frames for each counterfactual plot
# note that cmpgnMaxDisbursDvsR_mn is set to -1SD / Mean / +1SD for plot
# all other predictors (save one) set to the mean value
# imagining all districts somehow had a Democratic party incumbent each time
newDat_disbMinus1 <-
tibble(district2 = houseElectionDat_std_agg$district2,
nObs = rep(3, 435),
incumbentD_mn = houseElectionDat_std_agg$incumbentD_mn,
incumbentR_mn = houseElectionDat_std_agg$incumbentR_mn,
primaryUnopD_mn =
rep(mean(houseElectionDat_std_agg$primaryUnopD_mn), 435),
primaryUnopR_mn =
rep(mean(houseElectionDat_std_agg$primaryUnopR_mn), 435),
cmpgnMaxDisbursDvsR_mn =
rep(-1, 435),
primaryTotDvsR_mn =
rep(mean(houseElectionDat_std_agg$primaryTotDvsR_mn), 435),
log_popDensPerSqMi_mn =
rep(mean(houseElectionDat_std_agg$log_popDensPerSqMi_mn), 435),
medianAgeYr_mn =
rep(mean(houseElectionDat_std_agg$medianAgeYr_mn), 435),
pctEduGradProAge25plus_mn =
rep(mean(houseElectionDat_std_agg$pctEduGradProAge25plus_mn), 435),
pctInLaborForceUnemp_mn =
rep(mean(houseElectionDat_std_agg$pctInLaborForceUnemp_mn), 435),
perCapitaIncome_mn =
rep(mean(houseElectionDat_std_agg$perCapitaIncome_mn), 435),
pctWhiteAlone_mn =
rep(mean(houseElectionDat_std_agg$pctWhiteAlone_mn), 435) )
# imagining all districts somehow had a Republican party incumbent each time
newDat_disbMean <-
tibble(district2 = houseElectionDat_std_agg$district2,
nObs = rep(3, 435),
incumbentD_mn = houseElectionDat_std_agg$incumbentD_mn,
incumbentR_mn = houseElectionDat_std_agg$incumbentR_mn,
primaryUnopD_mn =
rep(mean(houseElectionDat_std_agg$primaryUnopD_mn), 435),
primaryUnopR_mn =
rep(mean(houseElectionDat_std_agg$primaryUnopR_mn), 435),
cmpgnMaxDisbursDvsR_mn =
rep(mean(houseElectionDat_std_agg$cmpgnMaxDisbursDvsR_mn), 435),
primaryTotDvsR_mn =
rep(mean(houseElectionDat_std_agg$primaryTotDvsR_mn), 435),
log_popDensPerSqMi_mn =
rep(mean(houseElectionDat_std_agg$log_popDensPerSqMi_mn), 435),
medianAgeYr_mn =
rep(mean(houseElectionDat_std_agg$medianAgeYr_mn), 435),
pctEduGradProAge25plus_mn =
rep(mean(houseElectionDat_std_agg$pctEduGradProAge25plus_mn), 435),
pctInLaborForceUnemp_mn =
rep(mean(houseElectionDat_std_agg$pctInLaborForceUnemp_mn), 435),
perCapitaIncome_mn =
rep(mean(houseElectionDat_std_agg$perCapitaIncome_mn), 435),
pctWhiteAlone_mn =
rep(mean(houseElectionDat_std_agg$pctWhiteAlone_mn), 435) )
# imagining all districts somehow had neither a Democratic party incumbent
# nor a Republican party incumbent each time
newDat_disbPlus1 <-
tibble(district2 = houseElectionDat_std_agg$district2,
nObs = rep(3, 435),
incumbentD_mn = houseElectionDat_std_agg$incumbentD_mn,
incumbentR_mn = houseElectionDat_std_agg$incumbentR_mn,
primaryUnopD_mn =
rep(mean(houseElectionDat_std_agg$primaryUnopD_mn), 435),
primaryUnopR_mn =
rep(mean(houseElectionDat_std_agg$primaryUnopR_mn), 435),
cmpgnMaxDisbursDvsR_mn =
rep(1, 435),
primaryTotDvsR_mn =
rep(mean(houseElectionDat_std_agg$primaryTotDvsR_mn), 435),
log_popDensPerSqMi_mn =
rep(mean(houseElectionDat_std_agg$log_popDensPerSqMi_mn), 435),
medianAgeYr_mn =
rep(mean(houseElectionDat_std_agg$medianAgeYr_mn), 435),
pctEduGradProAge25plus_mn =
rep(mean(houseElectionDat_std_agg$pctEduGradProAge25plus_mn), 435),
pctInLaborForceUnemp_mn =
rep(mean(houseElectionDat_std_agg$pctInLaborForceUnemp_mn), 435),
perCapitaIncome_mn =
rep(mean(houseElectionDat_std_agg$perCapitaIncome_mn), 435),
pctWhiteAlone_mn =
rep(mean(houseElectionDat_std_agg$pctWhiteAlone_mn), 435) )
predCount_disbMinus1 <- posterior_predict(m2.allFull_vIvS2_int1,
newdata = newDat_disbMinus1)
predCount_disbMean <- posterior_predict(m2.allFull_vIvS2_int1,
newdata = newDat_disbMean)
predCount_disbPlus1 <- posterior_predict(m2.allFull_vIvS2_int1,
newdata = newDat_disbPlus1)
# need to take the mean of each column for estimate, then SD
predCount_Dat <-
tibble(district2 = houseElectionDat_std_agg$district2,
generalWinD = houseElectionDat_std_agg$generalWinD,
incumbentD = houseElectionDat_std_agg$incumbentD_mn,
incumbentR = houseElectionDat_std_agg$incumbentR_mn,
predCt_disbMinus1SD = apply(predCount_disbMinus1, 2, mean),
predSD_disbMinus1SD = apply(predCount_disbMinus1, 2, sd),
predCt_disbMean = apply(predCount_disbMean, 2, mean),
predSD_disbMean = apply(predCount_disbMean, 2, sd),
predCt_disbPlus1SD = apply(predCount_disbPlus1, 2, mean),
predSD_disbPlus1SD = apply(predCount_disbPlus1, 2, sd) )
# see SO user hadley's (accepted) answer
# https://stackoverflow.com/questions/25925556/gather-multiple-sets-of-columns
predCount_Dat <-
predCount_Dat %>%
gather(key, value, -district2, -generalWinD, -incumbentR, -incumbentD) %>%
extract(key, into = c("metric", "cmpgnMaxDisbursDvsR"),
# ^ = start of string,
# ([:alpha:]*) = first group; match upper/lower alphabetic characters
# match to end of group
# \\_ = literal underscore; separates first and second group here
# ([a-zA-Z0-9]*) = second group; match upper/lower alphanumeric characters
# match to end of group
regex = "^([:alpha:]*)\\_([a-zA-Z0-9]*)") %>%
spread(metric, value) %>%
mutate(cmpgnMaxDisbursDvsR = factor(cmpgnMaxDisbursDvsR,
levels = c("disbMinus1SD",
"disbMean",
"disbPlus1SD")),
# maximum of predCt - 2 SDs and 0
predMin = if_else(predCt - 2*predSD < 0, 0, predCt - 2*predSD),
# minimum of predCt + 2 SDs and 3
predMax = if_else(predCt + 2*predSD > 3, 3, predCt + 2*predSD),
incFactor = factor(
if_else(incumbentR*3 == 4, "4R",
if_else(incumbentR*3 == 3 & incumbentD*3 == 0, "3R",
if_else(incumbentR*3 == 3 & incumbentD*3 == 1, "3R 1D",
if_else(incumbentR*3 == 2 & incumbentD*3 == 0, "2R",
if_else(incumbentR*3 == 2 & incumbentD*3 == 1, "2R 1D",
if_else(incumbentR*3 == 1 & incumbentD*3 == 0, "1R",
if_else(incumbentR*3 == 0 & incumbentD*3 == 0, "None",
if_else(incumbentR*3 == 1 & incumbentD*3 == 1, "1D 1R",
if_else(incumbentR*3 == 2 & incumbentD*3 == 2, "2R 2D",
if_else(incumbentR*3 == 0 & incumbentD*3 == 1, "1D",
if_else(incumbentR*3 == 1 & incumbentD*3 == 2, "2D 1R",
if_else(incumbentR*3 == 0 & incumbentD*3 == 2, "2D",
if_else(incumbentR*3 == 1 & incumbentD*3 == 3, "3D 1R",
if_else(incumbentR*3 == 0 & incumbentD*3 == 3, "3D",
if_else(incumbentD*3 == 4, "4D", "CHECK"))))))))))))))),
levels = c("4R", "3R", "3R 1D", "2R", "2R 1D", "1R",
"None", "1D 1R", "2R 2D", "1D", "2D 1R",
"2D", "3D 1R", "3D", "4D", "CHECK"),
ordered = T) )
library(RColorBrewer)
# take the first -6- reds; last two are too light to be distinguishing
plotReds <- rev(brewer.pal(8, name = "Reds")) # original order is light to dark
# take the last two values, and use a grey for "none"
plotPurps <- brewer.pal(4, name = "Purples")
# opposite of plotReds, take the last -6- blues
plotBlues <- brewer.pal(8, name = "Blues")
# note: plot will use all 6 Reds, no Grey (for "None"),
# 2 Purples (for equal # D/R), and all 5 of 6 Blues (no 4th / #6)
plotCols <- c(plotReds[1:6], plotPurps[3:4], plotBlues[c(3:5,7,8)])
predCount_Dat %>%
ggplot(aes(x = district2, y = predCt, ymin = predMin, ymax = predMax,
color = incFactor)) +
geom_pointrange() +
facet_wrap(generalWinD ~ cmpgnMaxDisbursDvsR, ncol = 3) +
scale_color_manual(values = plotCols) +
#scale_color_viridis_d() +
labs(title = "Counterfactual plot considering 'campaignMaxDisbursDvsR'; Pt Est \u00B1 2 SDs ",
subtitle = "Numbers: Actual # of Democratic candidates elected across the three elections \n Groups: -1SD DvsR spend level / Predictor mean / +1SD assigned for all 3 elections \n Colors: # incumbents by party across the three elections",
x = NULL, y = "Predicted # Democratic candidates elected", color = NULL,
caption = "Due to redistricting in 2012, some districts had >3 incumbents") +
theme_ds1() +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
This counterfactual plot suggests that, based on the model, spend differential (maximum reported $ spent for a Democratic campaign minus the maximum spent for a Republican campaign) makes a sizable difference in the predicted number of candidates elected for a party. This is especially pronounced/fluid for districts with less partisan lean in incumbents.
At least two things should be noted here:
Because the predictor,
cmpgnMaxDisbursDvsR_mn, is not centered at 0, the -1 SD (left) column implies a greater spend differential in favor of the Republican party candidate than the +1 SD (right) column does for a Democratic party candidate. If working at the sample-side level for the disaggregated “natural scale” data (houseElectionDatinstead ofhouseElectionDat_std_agg) is correct here, the -1SD column implies the Republican party candidate maximum outspends the Democratic party maximum by roughly $2 million ($2,027,666). Likewise, the +1SD column implies the Democratic party candidate maximum outspends the Republican party maximum by roughly $1.6 million ($1,584,997).“Massively outspend all the other party’s candidates” is not a tenable strategy, not only due to presumable reduced marginal return but also because on the Republican side this would be an aggregate differential of $882 million per House election period, and on the Democratic side this would be a spend differential of about $689.5 million. According to the
houseElectionDatdata set, the most spent by any single campaign in the three years considered was FL-18, where $230.5 million dollars were spent by campaigns alone (excluding PACs / Super PACs / etc.) in 2012. In other words, the differential in spending per two-year election cycle would need to be at least 3 times the total of the most profligate single campaign in the data set, in order to align with this counterfactual analysis and model assumptions.
Considering some district-level posteriors for a few districts
Let’s consider posterior distribution estimates for some district-level intercepts and slopes, and population-level standard deviations, and correlations.
For this I will consider the following districts:
- PA-7, the district (just barely) containing Villanova University and whose contorted profile for this period has been described as “Goofy kicking Donald Duck”9. There was a Republican party incumbent each election, and the incumbent won each time.
MA-7, the district containing approximately half of Boston and much of its surrounding suburbs. Somewhat the obverse of PA-7, Massachusetts’ 7th had a Democratic party incumbent each election, and this incumbent won each time.
NH-1, the district containing the southeastern portion of New Hampshire. This district went back-and-forth between parties, with the same Democratic candidate elected in 2012 and 2016, and a Republican candidate elected in 2014. There was a Republican incumbent in 2012 and (the same candidate) in 2016, while there was a Democratic incumbent in 2014.
A few population-level posteriors are also considered.
This first code chunk extracts some population level parameters and the district-level deviances from the population estimates for each district considered. It then computes the district-level value estimates for each parameter.
post_m2.allFull_v2_int1 <-
m2.allFull_vIvS2_int1 %>%
posterior_samples() %>%
# selecting only a subset of parameters
select(b_Intercept, b_incumbentD_mn, b_incumbentR_mn,
b_cmpgnMaxDisbursDvsR_mn,
sd_district2__Intercept,
sd_district2__incumbentD_mn, sd_district2__incumbentR_mn,
sd_district2__cmpgnMaxDisbursDvsR_mn,
cor_district2__Intercept__incumbentD_mn,
cor_district2__Intercept__incumbentR_mn,
cor_district2__Intercept__cmpgnMaxDisbursDvsR_mn,
cor_district2__incumbentD_mn__cmpgnMaxDisbursDvsR_mn,
cor_district2__incumbentR_mn__cmpgnMaxDisbursDvsR_mn,
"r_district2[MA-7,Intercept]", "r_district2[NH-1,Intercept]",
"r_district2[PA-7,Intercept]",
"r_district2[MA-7,incumbentD_mn]", "r_district2[NH-1,incumbentD_mn]",
"r_district2[PA-7,incumbentD_mn]",
"r_district2[MA-7,incumbentR_mn]", "r_district2[NH-1,incumbentR_mn]",
"r_district2[PA-7,incumbentR_mn]",
"r_district2[MA-7,cmpgnMaxDisbursDvsR_mn]",
"r_district2[NH-1,cmpgnMaxDisbursDvsR_mn]",
"r_district2[PA-7,cmpgnMaxDisbursDvsR_mn]" )
post_m2_v2_int1_sub <-
post_m2.allFull_v2_int1 %>%
transmute(`intercept_MA-7` = b_Intercept + `r_district2[MA-7,Intercept]`,
`incumbentD_MA-7` = b_incumbentD_mn +
`r_district2[MA-7,incumbentD_mn]`,
`incumbentR_MA-7` = b_incumbentR_mn +
`r_district2[MA-7,incumbentR_mn]`,
`cmpgnMaxDisbursDvsR_MA-7` =
b_cmpgnMaxDisbursDvsR_mn +
`r_district2[MA-7,cmpgnMaxDisbursDvsR_mn]`,
`intercept_NH-1` = b_Intercept + `r_district2[NH-1,Intercept]`,
`incumbentD_NH-1` = b_incumbentD_mn +
`r_district2[NH-1,incumbentD_mn]`,
`incumbentR_NH-1` = b_incumbentR_mn +
`r_district2[NH-1,incumbentR_mn]`,
`cmpgnMaxDisbursDvsR_NH-1` =
b_cmpgnMaxDisbursDvsR_mn +
`r_district2[NH-1,cmpgnMaxDisbursDvsR_mn]`,
`intercept_PA-7` = b_Intercept + `r_district2[PA-7,Intercept]`,
`incumbentD_PA-7` = b_incumbentD_mn +
`r_district2[PA-7,incumbentD_mn]`,
`incumbentR_PA-7` = b_incumbentR_mn +
`r_district2[PA-7,incumbentR_mn]`,
`cmpgnMaxDisbursDvsR_PA-7` =
b_cmpgnMaxDisbursDvsR_mn +
`r_district2[PA-7,cmpgnMaxDisbursDvsR_mn]`) %>%
# gather to narrow data format
gather() %>%
# extract variables based on district2
extract(key, into = c("par", "district2"),
# ^ = start of string,
# ([:alpha:]*) = first group; match upper/lower alphabetic characters
# match to end of group
# \\_ = literal underscore; separates first and second group here
# ([a-zA-Z0-9\\-]*) = second group; match upper/lower alphanumeric chr
# AND (escaped) dash symbol; match to end of group
regex = "^([:alpha:]*)\\_([a-zA-Z0-9\\-]*)") %>%
# there are 5,250 post-warmup MCMC iterations per parameter (12 pars here)
mutate(obs = rep(1:5250, 12) ) %>%
# spread the parameters into their own columns,
# now linked to which iter and district2; value is the stored value per entry
spread(par, value)Next we’ll plot each district’s estimates for b_Intercept (the baseline tendency to elect Democratic candidates) and each of b_incumbentD and b_incumbentR (the estimated marginal effect a given party having an incumbent in the district associated with electing a Democratic candidate in that district).
plot1 <-
post_m2_v2_int1_sub %>%
ggplot(aes(x = intercept, y = incumbentD, color = district2)) +
geom_point(alpha = 0.4) +
#geom_text(aes(label = district2, x = 0, y = 0)) +
labs(title = "Posterior estimates for 3 districts",
x = "Intercept", y = "incumbentD") +
scale_color_manual(values = c(color_set8[1:2], "orangered")) +
theme_ds1() + theme(legend.position = "top")
plot2 <-
post_m2.allFull_v2_int1 %>%
ggplot(aes(x = cor_district2__Intercept__incumbentD_mn)) +
geom_density(fill = "gray", alpha = 0.7) +
labs(title = "Posterior estimate for cor parameter",
subtitle = "District-level correlation?",
x = "cor_Intercept_incumbentD", y = "Density") +
theme_ds1()
plot3 <-
post_m2.allFull_v2_int1 %>%
ggplot(aes(x = sd_district2__Intercept)) +
geom_density(fill = "gray", alpha = 0.7) +
labs(title = "Posterior estimate for SD parameter",
subtitle = "Inter-district SD?",
x = "sd_Intercept", y = "Density") +
theme_ds1()
plot4 <-
post_m2.allFull_v2_int1 %>%
ggplot(aes(x = sd_district2__incumbentD_mn)) +
geom_density(fill = "gray", alpha = 0.7) +
labs(title = "Posterior estimate for SD parameter",
subtitle = "Inter-district SD?",
x = "sd_incumbentD", y = "Density") +
theme_ds1()
library(cowplot)
plot_grid(plot1, plot2, plot3, plot4, nrow = 2)
The quarter of plots above show three district-level bivariate posterior estimates (upper left) as well as three population-level posterior estimates for correlation between the two parameters (upper right) and standard deviations which I believe summarize district-level estimates of variance from the population-level estimate for each parameter (bottom).
Each of the three districts appear to have considerable uncertainty for the b_Intercept and b_incumbentD estimates, with the former centered over zero and the latter centered over 12 or so.
I think the correlation estimate’s being centered over zero may be related to brms’s using non-centered parametrization, but I don’t understand the concept well enough to be very confident in that.
For comparison, here is the same setup with incumbentR:
plot1 <-
post_m2_v2_int1_sub %>%
ggplot(aes(x = intercept, y = incumbentR, color = district2)) +
geom_point(alpha = 0.4) +
#geom_text(aes(label = district2, x = 0, y = 0)) +
labs(title = "Posterior estimates for 3 districts",
x = "Intercept", y = "incumbentR") +
scale_color_manual(values = c(color_set8[1:2], "orangered")) +
theme_ds1() + theme(legend.position = "top")
plot2 <-
post_m2.allFull_v2_int1 %>%
ggplot(aes(x = cor_district2__Intercept__incumbentR_mn)) +
geom_density(fill = "gray", alpha = 0.7) +
labs(title = "Posterior estimate for cor parameter",
subtitle = "District-level correlation?",
x = "cor_Intercept_incumbentR", y = "Density") +
theme_ds1()
plot3 <-
post_m2.allFull_v2_int1 %>%
ggplot(aes(x = sd_district2__Intercept)) +
geom_density(fill = "gray", alpha = 0.7) +
labs(title = "Posterior estimate for SD parameter",
subtitle = "Inter-district SD?",
x = "sd_Intercept", y = "Density") +
theme_ds1()
plot4 <-
post_m2.allFull_v2_int1 %>%
ggplot(aes(x = sd_district2__incumbentR_mn)) +
geom_density(fill = "gray", alpha = 0.7) +
labs(title = "Posterior estimate for SD parameter",
subtitle = "Inter-district SD?",
x = "sd_incumbentR", y = "Density") +
theme_ds1()
plot_grid(plot1, plot2, plot3, plot4, nrow = 2)
There appears to be roughly equal uncertainty for each district’s estimate of b_incumbentR in the opposite direction from b_incumbentD. The population-level sd_incumbentR estimate has more proportional density over larger values compared to sd_incumbentD, suggesting perhaps there is wider variability for the former among districts.
Finally, here are estimates for each of b_incumbentD and b_incumbentR with cmpgnMaxDisbursDvsR.
plot1 <-
post_m2_v2_int1_sub %>%
ggplot(aes(x = incumbentD, y = cmpgnMaxDisbursDvsR, color = district2)) +
geom_point(alpha = 0.4) +
#geom_text(aes(label = district2, x = 0, y = 0)) +
labs(title = "Posterior estimates for 3 districts",
x = "incumbentD", y = "cmpgnMaxDisbursDvsR") +
scale_color_manual(values = c(color_set8[1:2], "orangered")) +
theme_ds1() + theme(legend.position = "top")
plot2 <-
post_m2.allFull_v2_int1 %>%
ggplot(aes(x = cor_district2__incumbentD_mn__cmpgnMaxDisbursDvsR_mn)) +
geom_density(fill = "gray", alpha = 0.7) +
labs(title = "Posterior estimate for cor parameter",
subtitle = "District-level correlation?",
x = "cor_incumbentD_cmpgnMaxDisbursDvsR", y = "Density") +
theme_ds1()
plot3 <-
post_m2.allFull_v2_int1 %>%
ggplot(aes(x = sd_district2__incumbentD_mn)) +
geom_density(fill = "gray", alpha = 0.7) +
labs(title = "Posterior estimate for SD parameter",
subtitle = "Inter-district SD?",
x = "sd_incumbentD", y = "Density") +
theme_ds1()
plot4 <-
post_m2.allFull_v2_int1 %>%
ggplot(aes(x = sd_district2__cmpgnMaxDisbursDvsR_mn)) +
geom_density(fill = "gray", alpha = 0.7) +
labs(title = "Posterior estimate for SD parameter",
subtitle = "Inter-district SD?",
x = "sd_cmpgnMaxDisbursDvsR", y = "Density") +
theme_ds1()
plot_grid(plot1, plot2, plot3, plot4, nrow = 2)
plot1 <-
post_m2_v2_int1_sub %>%
ggplot(aes(x = incumbentR, y = cmpgnMaxDisbursDvsR, color = district2)) +
geom_point(alpha = 0.4) +
#geom_text(aes(label = district2, x = 0, y = 0)) +
labs(title = "Posterior estimates for 3 districts",
x = "incumbentR", y = "cmpgnMaxDisbursDvsR") +
scale_color_manual(values = c(color_set8[1:2], "orangered")) +
theme_ds1() + theme(legend.position = "top")
plot2 <-
post_m2.allFull_v2_int1 %>%
ggplot(aes(x = cor_district2__incumbentR_mn__cmpgnMaxDisbursDvsR_mn)) +
geom_density(fill = "gray", alpha = 0.7) +
labs(title = "Posterior estimate for cor parameter",
subtitle = "District-level correlation?",
x = "cor_incumbentR_cmpgnMaxDisbursDvsR", y = "Density") +
theme_ds1()
plot3 <-
post_m2.allFull_v2_int1 %>%
ggplot(aes(x = sd_district2__incumbentR_mn)) +
geom_density(fill = "gray", alpha = 0.7) +
labs(title = "Posterior estimate for SD parameter",
subtitle = "Inter-district SD?",
x = "sd_incumbentR", y = "Density") +
theme_ds1()
plot4 <-
post_m2.allFull_v2_int1 %>%
ggplot(aes(x = sd_district2__cmpgnMaxDisbursDvsR_mn)) +
geom_density(fill = "gray", alpha = 0.7) +
labs(title = "Posterior estimate for SD parameter",
subtitle = "Inter-district SD?",
x = "sd_cmpgnMaxDisbursDvsR", y = "Density") +
theme_ds1()
plot_grid(plot1, plot2, plot3, plot4, nrow = 2)
Example of Stan code for the interaction-including model
This is the unwieldy code that brms converts from the brms::brm() design model into Stan code. Makes me very, very glad that brms exists…
stancode(m2.allFull_vIvS2_int1)## // generated with brms 2.6.0
## functions {
## }
## data {
## int<lower=1> N; // total number of observations
## int Y[N]; // response variable
## int trials[N]; // number of trials
## int<lower=1> K; // number of population-level effects
## matrix[N, K] X; // population-level design matrix
## // data for group-level effects of ID 1
## int<lower=1> J_1[N];
## int<lower=1> N_1;
## int<lower=1> M_1;
## vector[N] Z_1_1;
## vector[N] Z_1_2;
## vector[N] Z_1_3;
## vector[N] Z_1_4;
## vector[N] Z_1_5;
## vector[N] Z_1_6;
## vector[N] Z_1_7;
## vector[N] Z_1_8;
## vector[N] Z_1_9;
## vector[N] Z_1_10;
## vector[N] Z_1_11;
## vector[N] Z_1_12;
## vector[N] Z_1_13;
## int<lower=1> NC_1;
## int prior_only; // should the likelihood be ignored?
## }
## transformed data {
## int Kc = K - 1;
## matrix[N, K - 1] Xc; // centered version of X
## vector[K - 1] means_X; // column means of X before centering
## for (i in 2:K) {
## means_X[i - 1] = mean(X[, i]);
## Xc[, i - 1] = X[, i] - means_X[i - 1];
## }
## }
## parameters {
## vector[Kc] b; // population-level effects
## real temp_Intercept; // temporary intercept
## vector<lower=0>[M_1] sd_1; // group-level standard deviations
## matrix[M_1, N_1] z_1; // unscaled group-level effects
## // cholesky factor of correlation matrix
## cholesky_factor_corr[M_1] L_1;
## }
## transformed parameters {
## // group-level effects
## matrix[N_1, M_1] r_1 = (diag_pre_multiply(sd_1, L_1) * z_1)';
## vector[N_1] r_1_1 = r_1[, 1];
## vector[N_1] r_1_2 = r_1[, 2];
## vector[N_1] r_1_3 = r_1[, 3];
## vector[N_1] r_1_4 = r_1[, 4];
## vector[N_1] r_1_5 = r_1[, 5];
## vector[N_1] r_1_6 = r_1[, 6];
## vector[N_1] r_1_7 = r_1[, 7];
## vector[N_1] r_1_8 = r_1[, 8];
## vector[N_1] r_1_9 = r_1[, 9];
## vector[N_1] r_1_10 = r_1[, 10];
## vector[N_1] r_1_11 = r_1[, 11];
## vector[N_1] r_1_12 = r_1[, 12];
## vector[N_1] r_1_13 = r_1[, 13];
## }
## model {
## vector[N] mu = temp_Intercept + Xc * b;
## for (n in 1:N) {
## mu[n] += r_1_1[J_1[n]] * Z_1_1[n] + r_1_2[J_1[n]] * Z_1_2[n] + r_1_3[J_1[n]] * Z_1_3[n] + r_1_4[J_1[n]] * Z_1_4[n] + r_1_5[J_1[n]] * Z_1_5[n] + r_1_6[J_1[n]] * Z_1_6[n] + r_1_7[J_1[n]] * Z_1_7[n] + r_1_8[J_1[n]] * Z_1_8[n] + r_1_9[J_1[n]] * Z_1_9[n] + r_1_10[J_1[n]] * Z_1_10[n] + r_1_11[J_1[n]] * Z_1_11[n] + r_1_12[J_1[n]] * Z_1_12[n] + r_1_13[J_1[n]] * Z_1_13[n];
## }
## // priors including all constants
## target += normal_lpdf(b | 0, 10);
## target += normal_lpdf(temp_Intercept | 0, 10);
## target += cauchy_lpdf(sd_1 | 0, 10)
## - 13 * cauchy_lccdf(0 | 0, 10);
## target += normal_lpdf(to_vector(z_1) | 0, 1);
## target += lkj_corr_cholesky_lpdf(L_1 | 1);
## // likelihood including all constants
## if (!prior_only) {
## target += binomial_logit_lpmf(Y | trials, mu);
## }
## }
## generated quantities {
## // actual population-level intercept
## real b_Intercept = temp_Intercept - dot_product(means_X, b);
## corr_matrix[M_1] Cor_1 = multiply_lower_tri_self_transpose(L_1);
## vector<lower=-1,upper=1>[NC_1] cor_1;
## // take only relevant parts of correlation matrix
## cor_1[1] = Cor_1[1,2];
## cor_1[2] = Cor_1[1,3];
## cor_1[3] = Cor_1[2,3];
## cor_1[4] = Cor_1[1,4];
## cor_1[5] = Cor_1[2,4];
## cor_1[6] = Cor_1[3,4];
## cor_1[7] = Cor_1[1,5];
## cor_1[8] = Cor_1[2,5];
## cor_1[9] = Cor_1[3,5];
## cor_1[10] = Cor_1[4,5];
## cor_1[11] = Cor_1[1,6];
## cor_1[12] = Cor_1[2,6];
## cor_1[13] = Cor_1[3,6];
## cor_1[14] = Cor_1[4,6];
## cor_1[15] = Cor_1[5,6];
## cor_1[16] = Cor_1[1,7];
## cor_1[17] = Cor_1[2,7];
## cor_1[18] = Cor_1[3,7];
## cor_1[19] = Cor_1[4,7];
## cor_1[20] = Cor_1[5,7];
## cor_1[21] = Cor_1[6,7];
## cor_1[22] = Cor_1[1,8];
## cor_1[23] = Cor_1[2,8];
## cor_1[24] = Cor_1[3,8];
## cor_1[25] = Cor_1[4,8];
## cor_1[26] = Cor_1[5,8];
## cor_1[27] = Cor_1[6,8];
## cor_1[28] = Cor_1[7,8];
## cor_1[29] = Cor_1[1,9];
## cor_1[30] = Cor_1[2,9];
## cor_1[31] = Cor_1[3,9];
## cor_1[32] = Cor_1[4,9];
## cor_1[33] = Cor_1[5,9];
## cor_1[34] = Cor_1[6,9];
## cor_1[35] = Cor_1[7,9];
## cor_1[36] = Cor_1[8,9];
## cor_1[37] = Cor_1[1,10];
## cor_1[38] = Cor_1[2,10];
## cor_1[39] = Cor_1[3,10];
## cor_1[40] = Cor_1[4,10];
## cor_1[41] = Cor_1[5,10];
## cor_1[42] = Cor_1[6,10];
## cor_1[43] = Cor_1[7,10];
## cor_1[44] = Cor_1[8,10];
## cor_1[45] = Cor_1[9,10];
## cor_1[46] = Cor_1[1,11];
## cor_1[47] = Cor_1[2,11];
## cor_1[48] = Cor_1[3,11];
## cor_1[49] = Cor_1[4,11];
## cor_1[50] = Cor_1[5,11];
## cor_1[51] = Cor_1[6,11];
## cor_1[52] = Cor_1[7,11];
## cor_1[53] = Cor_1[8,11];
## cor_1[54] = Cor_1[9,11];
## cor_1[55] = Cor_1[10,11];
## cor_1[56] = Cor_1[1,12];
## cor_1[57] = Cor_1[2,12];
## cor_1[58] = Cor_1[3,12];
## cor_1[59] = Cor_1[4,12];
## cor_1[60] = Cor_1[5,12];
## cor_1[61] = Cor_1[6,12];
## cor_1[62] = Cor_1[7,12];
## cor_1[63] = Cor_1[8,12];
## cor_1[64] = Cor_1[9,12];
## cor_1[65] = Cor_1[10,12];
## cor_1[66] = Cor_1[11,12];
## cor_1[67] = Cor_1[1,13];
## cor_1[68] = Cor_1[2,13];
## cor_1[69] = Cor_1[3,13];
## cor_1[70] = Cor_1[4,13];
## cor_1[71] = Cor_1[5,13];
## cor_1[72] = Cor_1[6,13];
## cor_1[73] = Cor_1[7,13];
## cor_1[74] = Cor_1[8,13];
## cor_1[75] = Cor_1[9,13];
## cor_1[76] = Cor_1[10,13];
## cor_1[77] = Cor_1[11,13];
## cor_1[78] = Cor_1[12,13];
## }Appendix 1: - Single-year model fitting and analysis
Model 1 - Bernoulli model, using 2012 data
As a first approximation / candidate model, I will use the following setup:
Outcome variable: generalWinD (D winner in general election for a given House district)
Model form: Bernoulli
Justification: Each district’s outcome is binary; either a Democratic party candidate won the seat, or they did not (in which case, for each of the three years considered, a Republican party candidate won).
Where possible, I will account for measurement uncertainty in the predictor variables - this applies to variables from the U.S. Census Bureau.
I will also make use of partial pooling with varying intercepts by district.
Note: For each of the following models, I made sure the population-level parameters all had Rhat values of no more than 1.02 (all but one had values of 1.01 or less, with a strong majority having Rhat values of 1.00 or less), and I also made sure there was good stationarity and mixing of the MCMC chains.
These reviews were conducted in the console using summary() and plot(); more comprehensive analysis follows the model-fitting code chunk for selecting the best-performing model(s) to proceed with.
Planned next steps: evaluate model fits by an Information Criterion, and relative model weights by the same, then evaluate model predictions on the fit data, and finally create and assess predictions for the next election’s data; lather, rinse, repeat
Fitting Model 1 models
The next code chunk is pretty lengthy, but it details each model’s components and their priors therein.
# intercept only model
m1.12_int <-
brm(generalWinD ~ 1,
data = houseElectionDat12_std, family = bernoulli(),
prior = c(set_prior("normal(0,3)", class = "Intercept") ),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.95),
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.12_int")
# model including incumbency predictors ONLY
m1.12_incmbt <-
brm(generalWinD ~ 1 + incumbentD + incumbentR,
data = houseElectionDat12_std, family = bernoulli(),
prior = c(set_prior("normal(0,3)", class = "Intercept"),
set_prior("normal(0,3)", class = "b", coef = "incumbentD"),
set_prior("normal(0,3)", class = "b", coef = "incumbentR") ),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99),
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.12_incmbt")
# model including primary unopposed predictors ONLY
m1.12_prmryUnop <-
brm(generalWinD ~ 1 + primaryUnopD + primaryUnopR,
data = houseElectionDat12_std, family = bernoulli(),
prior = c(set_prior("normal(0,3)", class = "Intercept"),
set_prior("normal(0,2)", class = "b") ),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99),
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.12_prmryUnop")
# model including difference in max campaign spend ONLY
m1.12_spnd <-
brm(generalWinD ~ 1 + cmpgnMaxDisbursDvsR,
data = houseElectionDat12_std, family = bernoulli(),
prior = c(set_prior("normal(0,3)", class = "Intercept"),
set_prior("normal(0,2)", class = "b") ),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.95),
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.12_spnd")
# model including difference in primary vote totals ONLY
m1.12_prmryTot <-
brm(generalWinD ~ 1 + primaryTotDvsR,
data = houseElectionDat12_std, family = bernoulli(),
prior = c(set_prior("normal(0,3)", class = "Intercept"),
set_prior("normal(0,2)", class = "b") ),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.95),
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.12_prmryTot")
# model including population density ONLY
m1.12_logPopDens <-
brm(generalWinD ~ 1 + log_popDensPerSqMi,
data = houseElectionDat12_std, family = bernoulli(),
prior = c(set_prior("normal(0,3)", class = "Intercept"),
set_prior("normal(0,2)", class = "b") ),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.95),
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.12_logPopDens")
# model including median age ONLY
m1.12_mdnAge <-
brm(generalWinD ~ 1 + me(medianAgeYr, medianAgeYr_SE),
data = houseElectionDat12_std, family = bernoulli(),
prior = c(set_prior("normal(0,3)", class = "Intercept"),
set_prior("normal(0,2)", class = "b"),
set_prior("normal(0,10)", class = "meanme"),
set_prior("student_t(3, 0, 10)", class = "sdme")),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99),
thin = 2,
save_mevars = TRUE,
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.12_mdnAge")
# model including education level ONLY
m1.12_hiEdu <-
brm(generalWinD ~ 1 + me(pctEduGradProAge25plus, pctEduGradProAge25plus_SE),
data = houseElectionDat12_std, family = bernoulli(),
prior = c(set_prior("normal(0,3)", class = "Intercept"),
set_prior("normal(0,2)", class = "b"),
set_prior("normal(0,10)", class = "meanme"),
set_prior("student_t(3, 0, 10)", class = "sdme") ),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99),
thin = 2,
save_mevars = TRUE,
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.12_hiEdu")
# model including unemployment and per capita income ONLY
m1.12_unempInc <-
brm(generalWinD ~ 1 + me(pctInLaborForceUnemp, pctInLaborForceUnemp_SE) +
me(perCapitaIncome, perCapitaIncome_SE),
data = houseElectionDat12_std, family = bernoulli(),
prior = c(set_prior("normal(0,3)", class = "Intercept"),
set_prior("normal(0,2)", class = "b"),
set_prior("normal(0,10)", class = "meanme"),
set_prior("student_t(3, 0, 10)", class = "sdme") ),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99, max_treedepth = 12),
thin = 2,
save_mevars = TRUE,
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.12_unempInc")
# model including % white alone ONLY
m1.12_wht <-
brm(generalWinD ~ 1 + me(pctWhiteAlone, pctWhiteAlone_SE),
data = houseElectionDat12_std, family = bernoulli(),
prior = c(set_prior("normal(0,3)", class = "Intercept"),
set_prior("normal(0,2)", class = "b"),
set_prior("normal(0,10)", class = "meanme"),
set_prior("student_t(3, 0, 10)", class = "sdme")),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99),
save_mevars = TRUE,
thin = 2,
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.12_wht")
# model including election/campaign predictors ONLY
m1.12_elecCmpgn <-
brm(generalWinD ~ 1 + incumbentD + incumbentR + primaryUnopD + primaryUnopR +
cmpgnMaxDisbursDvsR + primaryTotDvsR,
data = houseElectionDat12_std, family = bernoulli(),
prior = c(set_prior("normal(0,3)", class = "Intercept"),
set_prior("normal(0,2)", class = "b") ),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99),
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.12_elecCmpgn")
# model including demographic predictors ONLY
# no measurement error inclusion due to issue with fitting information criterion
m1.12_demo_noME <-
brm(generalWinD ~ 1 + log_popDensPerSqMi + medianAgeYr + pctEduGradProAge25plus +
pctInLaborForceUnemp + perCapitaIncome + pctWhiteAlone,
data = houseElectionDat12_std, family = bernoulli(),
prior = c(set_prior("normal(0,3)", class = "Intercept"),
set_prior("normal(0,2)", class = "b") ),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99, max_treedepth = 12),
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.12_demo_noME")
# model including ALL predictors
# no measurement error inclusion due to issue with fitting information criterion
m1.12_fullHouse_noME <-
brm(generalWinD ~ 1 + incumbentD + incumbentR + primaryUnopD + primaryUnopR +
cmpgnMaxDisbursDvsR + primaryTotDvsR + log_popDensPerSqMi
+ medianAgeYr + pctEduGradProAge25plus +
pctInLaborForceUnemp + perCapitaIncome + pctWhiteAlone,
data = houseElectionDat12_std, family = bernoulli(),
prior = c(set_prior("normal(0,3)", class = "Intercept"),
set_prior("normal(0,2)", class = "b") ),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99, max_treedepth = 12),
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.12_fullHouse_noME")Evaluating Model 1 fits
Here, I’ll consider each model’s Information Criterion estimate, and conduct a more in-depth analysis on the model(s) that emerge as the most promising. This analysis will include residuals plots (district-level error in prediction), counterfactual plots (model-implied prediction for theoretical values of one predictor, holding other values constant at some level), and combination plots (to confirm the model(s)’s Markov chains indicate a good model fit).
Provided these steps turn out well, I will then investigate using the model(s) to predict 2014 data, and then move on to investigating whether fitting the same set of models to 2014 data exhibits similar behavior and also performance on predicting 2016 data.
Looking into fitting m1.12_unempInc
Evaluating Model 1 fits - Information Criterion & Weights
Comparing models by WAIC to estimate out-of-sample prediction evaluation
While I am able to use Leave-one-out cross-validation with Pareto-smoothed importance sampling (PSIS-LOO) to estimate each model’s out-of-sample prediction performance, there is a programmatic error when applying PSIS-LOO for computing model weights via brms::model_weights(); therefore, I’m using WAIC as an approximation for model weighting instead.
Update 12/23/18: Just noting for potential later reference - with m1.12_spnd, observations 105 and 306 were considered problematic/influential in the context of PSIS-LOO. Ditto for m1.12_elecCmpgn and m1.12_fullHouse_noME, the other two models containing the variable for differential in maximum campaign reported spend for D-vs.-R by district. Those two observations are for FL-18 and OH-8, respectively (each for 2012). FL-18 has a cmpgnMaxDisbursDvsR value of roughly -$1.4mn, and OH-8 has a value of about $-2.1mn. The OH-8 value is the largest-margin value favoring Republican spend (and has the greatest absolute value) in the data set, and the FL-18 value is third (and has an absolute value also leading any Democrat-favoring spend value).
psisLOOCV_m1.12 <- loo(m1.12_int,
m1.12_incmbt, m1.12_prmryUnop, m1.12_spnd, m1.12_prmryTot,
m1.12_logPopDens, m1.12_mdnAge, m1.12_hiEdu, m1.12_unempInc,
m1.12_wht, m1.12_elecCmpgn, m1.12_demo_noME,
m1.12_fullHouse_noME,
compare = T, reloo = T)
modWts_m1.12 <-
model_weights(m1.12_int,
m1.12_incmbt, m1.12_prmryUnop, m1.12_spnd, m1.12_prmryTot,
m1.12_logPopDens, m1.12_mdnAge, m1.12_hiEdu,
m1.12_unempInc,
m1.12_wht, m1.12_elecCmpgn, m1.12_demo_noME,
m1.12_fullHouse_noME,
weights = "waic")Quick overview of the exhaustive between-models PSIS-LOO comparison:
psisLOOCV_m1.12## LOOIC SE
## m1.12_int 602.42 3.16
## m1.12_incmbt 244.01 27.63
## m1.12_prmryUnop 606.27 3.52
## m1.12_spnd 271.68 82.05
## m1.12_prmryTot 348.75 25.61
## m1.12_logPopDens 495.16 19.65
## m1.12_mdnAge 577.85 10.48
## m1.12_hiEdu 585.94 9.12
## m1.12_unempInc 490.84 20.01
## m1.12_wht 476.79 18.75
## m1.12_elecCmpgn 185.49 46.34
## m1.12_demo_noME 417.06 25.36
## m1.12_fullHouse_noME 181.55 42.00
## m1.12_int - m1.12_incmbt 358.41 27.58
## m1.12_int - m1.12_prmryUnop -3.85 1.33
## m1.12_int - m1.12_spnd 330.74 81.92
## m1.12_int - m1.12_prmryTot 253.67 25.77
## m1.12_int - m1.12_logPopDens 107.26 19.35
## m1.12_int - m1.12_mdnAge 24.57 10.02
## m1.12_int - m1.12_hiEdu 16.48 8.50
## m1.12_int - m1.12_unempInc 111.58 19.79
## m1.12_int - m1.12_wht 125.63 18.38
## m1.12_int - m1.12_elecCmpgn 416.93 46.30
## m1.12_int - m1.12_demo_noME 185.36 25.06
## m1.12_int - m1.12_fullHouse_noME 420.87 41.96
## m1.12_incmbt - m1.12_prmryUnop -362.26 27.48
## m1.12_incmbt - m1.12_spnd -27.67 78.24
## m1.12_incmbt - m1.12_prmryTot -104.74 35.08
## m1.12_incmbt - m1.12_logPopDens -251.15 31.94
## m1.12_incmbt - m1.12_mdnAge -333.84 28.80
## m1.12_incmbt - m1.12_hiEdu -341.93 28.83
## m1.12_incmbt - m1.12_unempInc -246.83 32.84
## m1.12_incmbt - m1.12_wht -232.78 31.33
## m1.12_incmbt - m1.12_elecCmpgn 58.52 38.87
## m1.12_incmbt - m1.12_demo_noME -173.05 34.78
## m1.12_incmbt - m1.12_fullHouse_noME 62.46 34.57
## m1.12_prmryUnop - m1.12_spnd 334.60 81.92
## m1.12_prmryUnop - m1.12_prmryTot 257.52 26.30
## m1.12_prmryUnop - m1.12_logPopDens 111.11 19.50
## m1.12_prmryUnop - m1.12_mdnAge 28.42 10.12
## m1.12_prmryUnop - m1.12_hiEdu 20.33 8.63
## m1.12_prmryUnop - m1.12_unempInc 115.43 19.87
## m1.12_prmryUnop - m1.12_wht 129.48 18.49
## m1.12_prmryUnop - m1.12_elecCmpgn 420.78 46.29
## m1.12_prmryUnop - m1.12_demo_noME 189.21 25.18
## m1.12_prmryUnop - m1.12_fullHouse_noME 424.72 41.96
## m1.12_spnd - m1.12_prmryTot -77.08 83.03
## m1.12_spnd - m1.12_logPopDens -223.48 83.56
## m1.12_spnd - m1.12_mdnAge -306.17 80.91
## m1.12_spnd - m1.12_hiEdu -314.27 82.23
## m1.12_spnd - m1.12_unempInc -219.17 84.57
## m1.12_spnd - m1.12_wht -205.11 82.76
## m1.12_spnd - m1.12_elecCmpgn 86.19 42.87
## m1.12_spnd - m1.12_demo_noME -145.38 84.01
## m1.12_spnd - m1.12_fullHouse_noME 90.13 48.30
## m1.12_prmryTot - m1.12_logPopDens -146.41 31.29
## m1.12_prmryTot - m1.12_mdnAge -229.10 27.29
## m1.12_prmryTot - m1.12_hiEdu -237.19 27.40
## m1.12_prmryTot - m1.12_unempInc -142.09 31.57
## m1.12_prmryTot - m1.12_wht -128.04 29.70
## m1.12_prmryTot - m1.12_elecCmpgn 163.27 46.62
## m1.12_prmryTot - m1.12_demo_noME -68.31 33.83
## m1.12_prmryTot - m1.12_fullHouse_noME 167.20 42.52
## m1.12_logPopDens - m1.12_mdnAge -82.69 20.78
## m1.12_logPopDens - m1.12_hiEdu -90.78 18.06
## m1.12_logPopDens - m1.12_unempInc 4.32 20.51
## m1.12_logPopDens - m1.12_wht 18.37 20.06
## m1.12_logPopDens - m1.12_elecCmpgn 309.67 49.56
## m1.12_logPopDens - m1.12_demo_noME 78.10 18.50
## m1.12_logPopDens - m1.12_fullHouse_noME 313.61 44.29
## m1.12_mdnAge - m1.12_hiEdu -8.09 14.70
## m1.12_mdnAge - m1.12_unempInc 87.01 21.14
## m1.12_mdnAge - m1.12_wht 101.06 16.21
## m1.12_mdnAge - m1.12_elecCmpgn 392.36 45.84
## m1.12_mdnAge - m1.12_demo_noME 160.79 23.69
## m1.12_mdnAge - m1.12_fullHouse_noME 396.30 41.26
## m1.12_hiEdu - m1.12_unempInc 95.10 20.70
## m1.12_hiEdu - m1.12_wht 109.15 20.45
## m1.12_hiEdu - m1.12_elecCmpgn 400.45 47.09
## m1.12_hiEdu - m1.12_demo_noME 168.88 24.47
## m1.12_hiEdu - m1.12_fullHouse_noME 404.39 42.56
## m1.12_unempInc - m1.12_wht 14.05 18.90
## m1.12_unempInc - m1.12_elecCmpgn 305.35 50.47
## m1.12_unempInc - m1.12_demo_noME 73.78 17.08
## m1.12_unempInc - m1.12_fullHouse_noME 309.29 44.86
## m1.12_wht - m1.12_elecCmpgn 291.30 48.32
## m1.12_wht - m1.12_demo_noME 59.73 16.62
## m1.12_wht - m1.12_fullHouse_noME 295.24 43.28
## m1.12_elecCmpgn - m1.12_demo_noME -231.57 50.89
## m1.12_elecCmpgn - m1.12_fullHouse_noME 3.94 9.63
## m1.12_demo_noME - m1.12_fullHouse_noME 235.51 44.75mods <- names(modWts_m1.12)
modWts_m1.12 %>%
as_tibble() %>%
mutate(model = mods) %>%
rename(WAIC_weight = value) %>%
select(model, WAIC_weight) %>%
arrange(desc(WAIC_weight)) %>%
mutate(WAIC_weight = formatC(WAIC_weight, format = "e", digits = 2)) %>%
kable(align = "c", booktabs = T, digits = 3) %>%
kable_styling(full_width = F)| model | WAIC_weight |
|---|---|
| m1.12_fullHouse_noME | 8.97e-01 |
| m1.12_elecCmpgn | 1.03e-01 |
| m1.12_incmbt | 9.15e-15 |
| m1.12_spnd | 1.33e-20 |
| m1.12_prmryTot | 1.63e-37 |
| m1.12_demo_noME | 2.43e-52 |
| m1.12_wht | 2.57e-65 |
| m1.12_unempInc | 2.35e-68 |
| m1.12_logPopDens | 2.64e-69 |
| m1.12_mdnAge | 2.92e-87 |
| m1.12_hiEdu | 5.11e-89 |
| m1.12_int | 1.35e-92 |
| m1.12_prmryUnop | 1.97e-93 |
While no model has exactly zero weight, the “full” model has the greatest share of model weight by such a margin that I will use it alone going forward.
Evaluating Model 1 fits - the “full” model
Let’s confirm the model fit well, and then evaluate prior estimates, and see if the full model could possibly be trimmed:
## NULLModel fits look good - each chain follows essentially the same density distribution for each parameter, and there is good mixing and stationarity for each chain.
Visualize posterior estimates:
It seems very plausible that pctWhiteAlone and perhaps perCapitaIncome do not provide any marginal predictive insight into the probability of a Democratic candidate winning House election over a Republican competitor.
Evaluating Model 1 fits - Improving the “full” model
After first fitting a model dropping pctWhiteAlone (m1.12_redHouse), then one dropping both pctWhiteAlone and perCapitaIncome (m1.12_redHouse2), and finally a model dropping these two predictors and also medianAgeYr (m1.12_redHouse3), it appears that this last model is the best candidate, with lower PSIS-LOOIC values and greater allocated model weight (this time using LOO was successful):
| model | weight |
|---|---|
| m1.12_redHouse3 | 0.585 |
| m1.12_redHouse2 | 0.372 |
| m1.12_redHouse | 0.026 |
| m1.12_fullHouse_noME | 0.017 |
Since the model excluding pctWhiteAlone, perCapitaIncome, and medianAgeYr performs so much better than the others, I’ll proceed with it for now.
Each posterior is considerably more concentrated than the respective posterior, with relatively slim plausibility of “zero” values in each case.
Evaluating Model 1 fits - In-sample performance
Let’s assess how well m1.12_redHouse3 predicts the winner for in-sample data.
# prediction of each district (columns) for each post-warmup iteration (rows)
postPred_m1.12_redHouse3 <- posterior_predict(m1.12_redHouse3)
# mean for each district (column of postPred_m1.12_redHouse3)
meanPred_m1.12_redHouse3 <- apply(postPred_m1.12_redHouse3, 2, mean)
# consider vote margin of victory
generalVoteDvsR <-
houseElectionDat %>%
filter(year == 2012) %>%
mutate(x = generalMaxD - generalMaxR) %>%
select(x) %>%
# standardize to use relative magnitude rather than raw counts
mutate(x = (x - mean(x)) / sd(x))
predVsObs12 <-
tibble(district2 = houseElectionDat12_std$district2,
generalWinD = houseElectionDat12_std$generalWinD,
meanPredWinD = meanPred_m1.12_redHouse3,
genlVoteDvsR_std = generalVoteDvsR$x)
nBins <- 100
cols <- c("darkred", "red",
"lightgrey",
"blue", "darkblue")
colGradient <- colorRampPalette(cols)
cut.cols <- colGradient(72)
cuts <- cut(predVsObs12$meanPredWinD, nBins)
names(cuts) <- sapply(cuts,function(t)
cut.cols[which(as.character(t) == levels(cuts))])
predVsObs12 %>%
ggplot(aes(x = meanPredWinD, fill = cut(meanPredWinD, nBins))) +
geom_histogram(bins = nBins) +
facet_wrap(. ~ generalWinD) +
labs(title = "Posterior prediction check for m1.12_redHouse3",
subtitle = "Faceted by Democratic candidate elected (0 / 1)",
x = "Mean proportion of model predictions, D candidate elected",
y = "Count") +
scale_fill_manual(values = cut.cols,
guide = F) +
theme_ds1()
rm(postPred_m1.12_redHouse3, generalVoteDvsR)Overall the model seems to perform pretty well on in-sample data, but there are some cases where a district was incorrectly estimated in 75% or more of the Markov Chain samples. Let’s check into these cases.
# compute proportion of cases correctly assigned to generalWinD = 0 / 1
# using <= 25% & 0 // >= 75% & 1 as cutoff for acceptable accuracy
predVsObs12 <-
predVsObs12 %>%
mutate(hitOrMiss = if_else((generalWinD == 0 & meanPredWinD <= 0.25) |
(generalWinD == 1 & meanPredWinD >= 0.75), "hit",
"miss"),
missMargin = abs(generalWinD - meanPredWinD))
# data show as counts
predVsObs12 %>%
group_by(generalWinD, hitOrMiss) %>%
summarize(nObs = n())## # A tibble: 4 x 3
## # Groups: generalWinD [?]
## generalWinD hitOrMiss nObs
## <int> <chr> <int>
## 1 0 hit 218
## 2 0 miss 16
## 3 1 hit 182
## 4 1 miss 19# for which districts did the model "miss" most?
predVsObs12 %>%
filter(hitOrMiss == "miss") %>%
arrange(desc(missMargin)) %>%
head(20) %>%
kable(align = "c") %>%
kable_styling(full_width = F, bootstrap_options = c("striped"))| district2 | generalWinD | meanPredWinD | genlVoteDvsR_std | hitOrMiss | missMargin |
|---|---|---|---|---|---|
| FL-18 | 1 | 0.0016 | -0.0002 | miss | 0.9984 |
| PA-12 | 0 | 0.9941 | -0.1283 | miss | 0.9941 |
| NY-27 | 0 | 0.9608 | -0.0649 | miss | 0.9608 |
| IL-17 | 1 | 0.0703 | 0.1591 | miss | 0.9297 |
| NC-8 | 0 | 0.9254 | -0.2389 | miss | 0.9254 |
| NH-1 | 1 | 0.0924 | 0.1037 | miss | 0.9076 |
| NY-24 | 1 | 0.0956 | 0.1319 | miss | 0.9044 |
| FL-2 | 0 | 0.8342 | -0.1889 | miss | 0.8342 |
| TX-23 | 1 | 0.1698 | 0.0675 | miss | 0.8302 |
| KY-6 | 0 | 0.8211 | -0.1285 | miss | 0.8211 |
| NY-18 | 1 | 0.1939 | 0.0832 | miss | 0.8061 |
| NH-2 | 1 | 0.2353 | 0.1347 | miss | 0.7647 |
| MN-8 | 1 | 0.2500 | 0.2768 | miss | 0.7500 |
| OK-2 | 0 | 0.7171 | -0.4645 | miss | 0.7171 |
| IL-10 | 1 | 0.3081 | 0.0131 | miss | 0.6919 |
| CA-36 | 1 | 0.3315 | 0.0967 | miss | 0.6685 |
| NC-7 | 1 | 0.3562 | -0.0119 | miss | 0.6438 |
| IL-11 | 1 | 0.3574 | 0.3905 | miss | 0.6426 |
| CA-21 | 0 | 0.6372 | -0.1872 | miss | 0.6372 |
| FL-26 | 1 | 0.3860 | 0.2339 | miss | 0.6140 |
The output above shows that the model was “fooled” by about 45 districts. My rule here is the model had a bad “miss” if the mean of a district’s iterations indicate 75% or more iterations predicting a Democratic candidate winner when a Republican candidate won, OR 25% or fewer iterations predicting a Democratic candidate winner when a Democratic candidate won.
Looking at things a little more closely, each of the 20 districts with the greatest “miss margin” (discrepancy between mean predicted outcome and the actual outcome) had a Democratic-versus-Republican vote margin within a half standard deviation of the sample mean. I may very well be misreading things, but that seems to be an indication of the model not so much falling flat on its face (metaphorically speaking) as simply missing some underlying signal that would aid prediction accuracy.
If time allows, I will dig deeper into this, and at the very least see if the model fit to all three years’ data fares any better.
Evaluating model performance
At the 12/6/18 check-in with Dr. Frey, he suggested one way of evaluating the model’s predictive performance as checking the relative proportion of candidates with a given range of mean predicted success/failure and seeing how that compares to reality.
For example, if the model has a given district with a mean proportion of about 0.7-0.8 of all simulations resulting in the Democratic candidate being elected, what proportion of the time is the model correct? How about for cases where the model “predicts” the Republican candidate will be elected in about 70-80% of simulations?
sum(predVsObs12$meanPredWinD >= 0.7 & predVsObs12$meanPredWinD <= 0.8 &
predVsObs12$generalWinD == 1) /
sum(predVsObs12$meanPredWinD >= 0.7 & predVsObs12$meanPredWinD <= 0.8)## [1] 0.9sum(predVsObs12$meanPredWinD >= 0.50 & predVsObs12$meanPredWinD <= 0.60 &
predVsObs12$generalWinD == 1) /
sum(predVsObs12$meanPredWinD >= 0.50 & predVsObs12$meanPredWinD <= 0.60)## [1] 0.5sum(predVsObs12$meanPredWinD >= 0.40 & predVsObs12$meanPredWinD <= 0.50 &
predVsObs12$generalWinD == 0) /
sum(predVsObs12$meanPredWinD >= 0.40 & predVsObs12$meanPredWinD <= 0.50)## [1] 0.4sum(predVsObs12$meanPredWinD >= 0.2 & predVsObs12$meanPredWinD <= 0.3 &
predVsObs12$generalWinD == 0) /
sum(predVsObs12$meanPredWinD >= 0.2 & predVsObs12$meanPredWinD <= 0.3)## [1] 0.8462The model seems to “get it right” roughly as often as one would hope, with better performance in cases with a more pronounced partisan lean.
Where the model predicts a Democratic candidate’s being elected in 70-80% of simulations, a Democrat “wins” about 85% of the time. For cases where the model predicts a Democrat will be successful in 50-60% of simulations, they are successful about 50% of the time. For model-predicted “40-50% Democratic likelihood” cases, a Republican candidate is elected in about 40% of cases. And in cases where a Democratic party candidate is elected in 20-30% of simulations, the Republican candidate is elected 75% of the time.
A histogram of the “miss margin” might be useful to see as well - though this is essentially the aggregated version of the prediction plot.
predVsObs12 %>%
ggplot(aes(x = missMargin)) +
geom_histogram(bins = 100, color = "white", fill = color_set8[6]) +
theme_ds1()
Digging deeper into model “misses”
Perhaps it will be useful to put up some quick plots of the model predictors, and potential predictors not currently in the model, to see where things go astray and where the model is working well.
Note: This section was more or less set aside - nothing really jumped out to me, and I will need to carry out Dr. Frey’s suggestion to consider relative proportion of the “bad-miss” cases vs. that of the “non-/slight-miss” cases for these predictors to really identify which predictors may enable better classification of the current “bad misses”.
houseElectionDat12_std <-
houseElectionDat12_std %>%
mutate(m1.12_redHouse3_missMargin = predVsObs12$missMargin)
houseElectionDat12_std %>%
ggplot(aes(x = cmpgnMaxDisbursDvsR, y = m1.12_redHouse3_missMargin)) +
geom_vline(xintercept = 0, color = "gray") +
geom_point(aes(color = m1.12_redHouse3_missMargin)) +
scale_color_viridis_c(direction = -1) +
theme_ds1()
Evaluating Model 1 fits - Out-of-sample performance (2014 election)
Now the rubber meets the road - let’s see how the model fares in predicting the following House election.
# estimate is the mean predicted response (generalWinD)
# for median, use robust = T
pred14From12 <- predict(m1.12_redHouse3, newdata = houseElectionDat14_std)
# evaluate with same general approach as 2012 in-sample posterior pred check
# consider vote margin of victory
generalVoteDvsR <-
houseElectionDat %>%
filter(year == 2014) %>%
mutate(x = generalMaxD - generalMaxR) %>%
select(x) %>%
# standardize to use relative magnitude rather than raw counts
mutate(x = (x - mean(x)) / sd(x))
predVsObs14From12 <-
tibble(district2 = houseElectionDat14_std$district2,
generalWinD = houseElectionDat14_std$generalWinD,
meanPredWinD = pred14From12[,"Estimate"],
genlVoteDvsR_std = generalVoteDvsR$x)
predVsObs14From12 %>%
ggplot(aes(x = meanPredWinD, fill = cut(meanPredWinD, nBins))) +
geom_histogram(bins = nBins) +
facet_wrap(. ~ generalWinD) +
labs(title = "Prediction check for m1.12_redHouse3 using 2014 data",
subtitle = "Faceted by Democratic candidate elected (0 / 1)",
x = "Mean proportion of model predictions, D candidate elected",
y = "Count",
caption = "Note: model fit to 2012 election data") +
scale_fill_manual(values = cut.cols,
guide = F) +
theme_ds1()
# compute proportion of cases correctly assigned to generalWinD = 0 / 1
# using <= 25% & 0 // >= 75% & 1 as cutoff for acceptable accuracy
predVsObs14From12 <-
predVsObs14From12 %>%
mutate(hitOrMiss = if_else((generalWinD == 0 & meanPredWinD <= 0.25) |
(generalWinD == 1 & meanPredWinD >= 0.75), "hit",
"miss"),
missMargin = abs(generalWinD - meanPredWinD))
# data show as counts
predVsObs14From12 %>%
group_by(generalWinD, hitOrMiss) %>%
summarize(nObs = n())## # A tibble: 4 x 3
## # Groups: generalWinD [?]
## generalWinD hitOrMiss nObs
## <int> <chr> <int>
## 1 0 hit 219
## 2 0 miss 28
## 3 1 hit 182
## 4 1 miss 6# for which districts did the model "miss" most?
predVsObs14From12 %>%
filter(hitOrMiss == "miss") %>%
arrange(desc(missMargin)) %>%
head(20) %>%
kable(align = "c") %>%
kable_styling(full_width = F, bootstrap_options = c("striped"))| district2 | generalWinD | meanPredWinD | genlVoteDvsR_std | hitOrMiss | missMargin |
|---|---|---|---|---|---|
| WV-3 | 0 | 0.9811 | -0.0557 | miss | 0.9811 |
| NE-2 | 1 | 0.0200 | 0.2188 | miss | 0.9800 |
| IL-10 | 0 | 0.9629 | 0.0789 | miss | 0.9629 |
| NY-24 | 0 | 0.9504 | -0.3621 | miss | 0.9504 |
| NV-4 | 0 | 0.9481 | 0.0952 | miss | 0.9481 |
| IL-12 | 0 | 0.9129 | -0.1504 | miss | 0.9129 |
| TX-23 | 0 | 0.9030 | 0.1111 | miss | 0.9030 |
| NY-1 | 0 | 0.8955 | -0.0595 | miss | 0.8955 |
| VA-10 | 0 | 0.8510 | -0.3328 | miss | 0.8510 |
| GA-12 | 0 | 0.8057 | -0.0667 | miss | 0.8057 |
| FL-26 | 0 | 0.7730 | 0.0806 | miss | 0.7730 |
| FL-2 | 1 | 0.2479 | 0.1807 | miss | 0.7521 |
| AZ-2 | 0 | 0.7507 | 0.1410 | miss | 0.7507 |
| KY-6 | 0 | 0.7461 | -0.5069 | miss | 0.7461 |
| PA-6 | 0 | 0.6887 | -0.2108 | miss | 0.6887 |
| MN-7 | 1 | 0.3546 | 0.4157 | miss | 0.6454 |
| MN-8 | 1 | 0.4327 | 0.1926 | miss | 0.5673 |
| NH-1 | 0 | 0.5343 | 0.0275 | miss | 0.5343 |
| NY-19 | 0 | 0.5119 | -0.6394 | miss | 0.5119 |
| KY-5 | 0 | 0.4880 | -1.4946 | miss | 0.4880 |
The 2012-fit reduced model seems disproportionately less accurate in predicting Republican candidate wins for 2014 data. As stated previously, I am defining a “miss” as a mean predicted win outcome giving the Republican candidate \(\le\) 25% likelihood of success when a Republican is actually elected, or a mean predicted outcome giving the Democratic candidate \(\le\) 25% likelihood of success when a Democrat is actually elected.
By such a metric, the model’s predictions “miss” in about 11% of all cases where a Republican candidate won, versus about 5% of cases where a Democratic candidate won.
Model 1 - Bernoulli model, using 2014 data
Now let’s see how similar posterior distribution estimates are when fitting the “full” model to 2014 data, and if the 2014-fit “full” model can be similarly pared down as the 2012 model was.
m1.14_fullHouse_noME <-
brm(generalWinD ~ 1 + incumbentD + incumbentR + primaryUnopD + primaryUnopR +
cmpgnMaxDisbursDvsR + primaryTotDvsR + log_popDensPerSqMi
+ medianAgeYr + pctEduGradProAge25plus +
pctInLaborForceUnemp + perCapitaIncome + pctWhiteAlone,
data = houseElectionDat14_std, family = bernoulli(),
prior = c(set_prior("normal(0,3)", class = "Intercept"),
set_prior("normal(0,2)", class = "b") ),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99, max_treedepth = 12),
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.14_fullHouse_noME")
Interestingly, while the 2012-fit “full” model had pctWhiteAlone, perCapitaIncome, and medianAgeYr as predictors whose posterior distribution gave serious plausibility to zero-effect estimates, for the 2014-fit “full” model it is pctEduGradProAge25plus, primaryUnopR, and primaryUnopD that have posteriors suggesting their removal might not adversely impact the model’s in-sample predictive capabilities.
Just for so, let’s see how the 2014-fit model performs when predicting on in-sample data with all variables versus when those three variables are removed.
# prediction of 2014 model in-sample, "full" model
# prediction of each district (columns) for each post-warmup iteration (rows)
postPred_m1.14_fullHouse <- posterior_predict(m1.14_fullHouse_noME)
# mean for each district (column of postPred_m1.14_redHouse)
meanPred_m1.14_fullHouse <- apply(postPred_m1.14_fullHouse, 2, mean)
# reduced version of "full" 2014-fit model
m1.14_redHouse <-
update(m1.14_fullHouse_noME,
. ~ . -primaryUnopD -primaryUnopR -pctEduGradProAge25plus)## Start sampling# prediction of 2014 model in-sample, "reduced" model
# prediction of each district (columns) for each post-warmup iteration (rows)
postPred_m1.14_redHouse <- posterior_predict(m1.14_redHouse)
# mean for each district (column of postPred_m1.12_redHouse3)
meanPred_m1.14_redHouse <- apply(postPred_m1.14_redHouse, 2, mean)
# consider vote margin of victory
generalVoteDvsR <-
houseElectionDat %>%
filter(year == 2014) %>%
mutate(x = generalMaxD - generalMaxR) %>%
select(x) %>%
# standardize to use relative magnitude rather than raw counts
mutate(x = (x - mean(x)) / sd(x))
predVsObs14 <-
tibble(district2 = houseElectionDat14_std$district2,
generalWinD = houseElectionDat14_std$generalWinD,
meanPredWinDFull = meanPred_m1.14_fullHouse,
meanPredWinDRed = meanPred_m1.14_redHouse,
genlVoteDvsR_std = generalVoteDvsR$x)
predVsObs14 %>%
ggplot(aes(x = meanPredWinDFull, fill = cut(meanPredWinDFull, nBins))) +
geom_histogram(bins = nBins) +
facet_wrap(. ~ generalWinD) +
labs(title = "Posterior prediction check for m1.14_fullHouse_noME",
subtitle = "Faceted by Democratic candidate elected (0 / 1)",
x = "Mean proportion of model predictions, D candidate elected",
y = "Count") +
scale_fill_manual(values = cut.cols,
guide = F) +
theme_ds1()
predVsObs14 %>%
ggplot(aes(x = meanPredWinDRed, fill = cut(meanPredWinDRed, nBins))) +
geom_histogram(bins = nBins) +
facet_wrap(. ~ generalWinD) +
labs(title = "Posterior prediction check for m1.14_redHouse",
subtitle = "Faceted by Democratic candidate elected (0 / 1)",
x = "Mean proportion of model predictions, D candidate elected",
y = "Count") +
scale_fill_manual(values = cut.cols,
guide = F) +
theme_ds1()
# comparison of the two models' predictions
predVsObs14 %>%
ggplot() +
geom_density(aes(x = meanPredWinDFull), fill = color_set8[1], alpha = 0.4,
adjust = 0.2) +
geom_density(aes(x = meanPredWinDRed), fill = color_set8[2], alpha = 0.4,
adjust = 0.2) +
labs(title = "Comparing posterior prediction densities, 2014 models",
subtitle = "Full (dk blue) vs Reduced (lt blue)",
x = "Mean proportion of model predictions, D candidate elected",
y = "Density") +
theme_ds1()
The two models are impressively similar, and yield largely accurate inferences on election outcomes, as did the 2012-fit model. Fit diagnostics (trace plots, chain density estimates, and R-hat values) were good for both the full and reduced models again.
Out-of-sample performance - 2016 election
As with the 2012 model, let’s see how well the 2014 (reduced) model fares in predicting 2016 election outcomes.
# estimate is the mean predicted response (generalWinD)
# for median, use robust = T
pred16From14 <- predict(m1.14_redHouse, newdata = houseElectionDat16_std)
# evaluate with same general approach as 2012 in-sample posterior pred check
# consider vote margin of victory
generalVoteDvsR <-
houseElectionDat %>%
filter(year == 2016) %>%
mutate(x = generalMaxD - generalMaxR) %>%
select(x) %>%
# standardize to use relative magnitude rather than raw counts
mutate(x = (x - mean(x)) / sd(x))
predVsObs16From14 <-
tibble(district2 = houseElectionDat16_std$district2,
generalWinD = houseElectionDat16_std$generalWinD,
meanPredWinD = pred16From14[,"Estimate"],
genlVoteDvsR_std = generalVoteDvsR$x)
predVsObs16From14 %>%
ggplot(aes(x = meanPredWinD, fill = cut(meanPredWinD, nBins))) +
geom_histogram(bins = nBins) +
facet_wrap(. ~ generalWinD) +
labs(title = "Prediction check for m1.14_redHouse using 2016 data",
subtitle = "Faceted by Democratic candidate elected (0 / 1)",
x = "Mean proportion of model predictions, D candidate elected",
y = "Count",
caption = "Note: model fit to 2014 election data") +
scale_fill_manual(values = cut.cols,
guide = F) +
theme_ds1()
# compute proportion of cases correctly assigned to generalWinD = 0 / 1
# using <= 25% & 0 // >= 75% & 1 as cutoff for acceptable accuracy
predVsObs16From14 <-
predVsObs16From14 %>%
mutate(hitOrMiss = if_else((generalWinD == 0 & meanPredWinD <= 0.25) |
(generalWinD == 1 & meanPredWinD >= 0.75), "hit",
"miss"),
missMargin = abs(generalWinD - meanPredWinD))
# data show as counts
predVsObs16From14 %>%
group_by(generalWinD, hitOrMiss) %>%
summarize(nObs = n())## # A tibble: 4 x 3
## # Groups: generalWinD [?]
## generalWinD hitOrMiss nObs
## <int> <chr> <int>
## 1 0 hit 236
## 2 0 miss 5
## 3 1 hit 174
## 4 1 miss 20# for which districts did the model "miss" most?
predVsObs16From14 %>%
filter(hitOrMiss == "miss") %>%
arrange(desc(missMargin)) %>%
head(20) %>%
kable(align = "c") %>%
kable_styling(full_width = F, bootstrap_options = c("striped"))| district2 | generalWinD | meanPredWinD | genlVoteDvsR_std | hitOrMiss | missMargin |
|---|---|---|---|---|---|
| NH-1 | 1 | 0.0043 | 0.0804 | miss | 0.9957 |
| FL-7 | 1 | 0.0073 | 0.1257 | miss | 0.9927 |
| FL-13 | 1 | 0.0089 | 0.1506 | miss | 0.9911 |
| NV-4 | 1 | 0.0210 | 0.1273 | miss | 0.9790 |
| NJ-5 | 1 | 0.0545 | 0.1615 | miss | 0.9455 |
| AZ-1 | 1 | 0.0735 | 0.2065 | miss | 0.9265 |
| NE-2 | 0 | 0.8733 | 0.0133 | miss | 0.8733 |
| MN-8 | 1 | 0.1662 | 0.0575 | miss | 0.8338 |
| NV-3 | 1 | 0.1667 | 0.0731 | miss | 0.8333 |
| FL-18 | 0 | 0.8008 | -0.2782 | miss | 0.8008 |
| IL-10 | 1 | 0.2169 | 0.1615 | miss | 0.7831 |
| CA-24 | 1 | 0.3333 | 0.2128 | miss | 0.6667 |
| FL-9 | 1 | 0.3611 | 0.4518 | miss | 0.6389 |
| FL-26 | 0 | 0.6219 | -0.2256 | miss | 0.6219 |
| MN-7 | 1 | 0.4090 | 0.1756 | miss | 0.5910 |
| MN-2 | 0 | 0.5855 | -0.0125 | miss | 0.5855 |
| VA-4 | 1 | 0.4243 | 0.4804 | miss | 0.5757 |
| ME-1 | 1 | 0.5262 | 0.5496 | miss | 0.4738 |
| IL-2 | 1 | 0.5832 | 1.4586 | miss | 0.4168 |
| TX-15 | 1 | 0.5838 | 0.3225 | miss | 0.4162 |
The 2014-fit reduced model is disproportionately less accurate in predicting Democratic candidate wins for 2016 data, which is opposite of how the 2012-fit model fared in predicting 2014 outcomes. Now the model “missed” in 14% of cases where a Democratic candidate was elected, versus about 4% of cases where a Republican candidate won.
Model 1 - Bernoulli model, using 2016 data
Now let’s see how similar posterior distribution estimates are when fitting the “full” model to 2016 data, and if the 2016-fit “full” model can be similarly pared down as the 2012 and 2014 models were.
Note: For greater comparability, posterior coefficient distribution estimates for each year’s full model will be compared in the next section.
m1.16_fullHouse_noME <-
brm(generalWinD ~ 1 + incumbentD + incumbentR + primaryUnopD + primaryUnopR +
cmpgnMaxDisbursDvsR + primaryTotDvsR + log_popDensPerSqMi
+ medianAgeYr + pctEduGradProAge25plus +
pctInLaborForceUnemp + perCapitaIncome + pctWhiteAlone,
data = houseElectionDat16_std, family = bernoulli(),
prior = c(set_prior("normal(0,3)", class = "Intercept"),
set_prior("normal(0,2)", class = "b") ),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99, max_treedepth = 12),
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.16_fullHouse_noME")
The 2016 “full” model once again does not directly accord with prior-year “full” models, though it seems slightly more aligned with its 2012 counterpart compared to 2014.
For the 2016 model, medianAgeYr seems to very plausibly not contribute any marginal predictive information when all other predictors are accounted for. perCapitaIncome could also potentially be removed at little cost to predictive accuracy. The intercept appears to have a posterior distribution estimate most notably similar to that of medianAgeYr - I believe that this implies the predictors in the model are contributing the lion’s share of the marginal predictive
Unlike the 2012-fit “full” model, now primaryUnopD appears to provide a healthy amount of marginal predictive utility.
As before, let’s see how the 2016-fit model performs when predicting on in-sample data with all variables versus when those three variables are removed.
# prediction of 2016 model in-sample, "full" model
# prediction of each district (columns) for each post-warmup iteration (rows)
postPred_m1.16_fullHouse <- posterior_predict(m1.16_fullHouse_noME)
# mean for each district (column of postPred_m1.16_redHouse)
meanPred_m1.16_fullHouse <- apply(postPred_m1.16_fullHouse, 2, mean)
# reduced version of "full" 2014-fit model
m1.16_redHouse <-
update(m1.16_fullHouse_noME,
. ~ . -medianAgeYr -perCapitaIncome)## Start sampling# prediction of 2016 model in-sample, "reduced" model
# prediction of each district (columns) for each post-warmup iteration (rows)
postPred_m1.16_redHouse <- posterior_predict(m1.16_redHouse)
# mean for each district (column of postPred_m1.16_redHouse)
meanPred_m1.16_redHouse <- apply(postPred_m1.16_redHouse, 2, mean)
# consider vote margin of victory
generalVoteDvsR <-
houseElectionDat %>%
filter(year == 2016) %>%
mutate(x = generalMaxD - generalMaxR) %>%
select(x) %>%
# standardize to use relative magnitude rather than raw counts
mutate(x = (x - mean(x)) / sd(x))
predVsObs16 <-
tibble(district2 = houseElectionDat16_std$district2,
generalWinD = houseElectionDat16_std$generalWinD,
meanPredWinDFull = meanPred_m1.16_fullHouse,
meanPredWinDRed = meanPred_m1.16_redHouse,
genlVoteDvsR_std = generalVoteDvsR$x)
predVsObs16 %>%
ggplot(aes(x = meanPredWinDFull, fill = cut(meanPredWinDFull, nBins))) +
geom_histogram(bins = nBins) +
facet_wrap(. ~ generalWinD) +
labs(title = "Posterior prediction check for m1.16_fullHouse_noME",
subtitle = "Faceted by Democratic candidate elected (0 / 1)",
x = "Mean proportion of model predictions, D candidate elected",
y = "Count") +
scale_fill_manual(values = cut.cols,
guide = F) +
theme_ds1()
predVsObs16 %>%
ggplot(aes(x = meanPredWinDRed, fill = cut(meanPredWinDRed, nBins))) +
geom_histogram(bins = nBins) +
facet_wrap(. ~ generalWinD) +
labs(title = "Posterior prediction check for m1.16_redHouse",
subtitle = "Faceted by Democratic candidate elected (0 / 1)",
x = "Mean proportion of model predictions, D candidate elected",
y = "Count") +
scale_fill_manual(values = cut.cols,
guide = F) +
theme_ds1()
# comparison of the two models' predictions
predVsObs16 %>%
ggplot() +
geom_density(aes(x = meanPredWinDFull), fill = color_set8[5], alpha = 0.4,
adjust = 0.2) +
geom_density(aes(x = meanPredWinDRed), fill = color_set8[6], alpha = 0.4,
adjust = 0.2) +
labs(title = "Comparing posterior prediction densities, 2016 models",
subtitle = "Full (tan) vs Reduced (yellow)",
x = "Mean proportion of model predictions, D candidate elected",
y = "Density") +
theme_ds1()
As with the 2014 models, the 2016 “full” and reduced models are very strongly similar, as would easily be expected.
Appendix 2 - Bernoulli outcome models - good use of partial pooling, or island of misfit toy models?
Note: After going through the full analysis of this section, I’ve concluded that the analysis is pretty much bunk, but I’ll keep in Appendix 2 for reference.
The following models are fit treating the outcome variable (0 / 1 for generalWinD) as a Bernoulli random variable - in other words, each district-level election has a single trial with some unknown probability of “success” (here, the election of a Democratic party candidate), to be estimated from the Bayesian regression.
After going back and forth on this (in consultation with chapter 12 of Statistical Rethinking, which discusses multilevel models), I’m proceeding with the Bernoulli-outcome models as the “good” models for deeper analysis and consideration in this project, rather than the next section’s Binomial-outcome models (which use the same data, just aggregated at the district level for the three elections).
I believe this section’s models are the better application of multilevel modeling for one distinct reason. Before getting to that, I believe that using district2 (the variable identifying unique districts, e.g. “PA-1” being Pennsylvania’s 1st district) as the “cluster-identifier” is correct in either the Bernoulli or the Binomial case, and either case’s varying-effects model will be able to partially pool varying intercept / varying slopes inferences among districts identified as having similar underlying behavior (association between predictor and outcome values).
What makes the Bernoulli case, explored in this section, a potentially better application of partial pooling is that, because each district-level election is its own observation (rather than 1/3 of an aggregated district-level observation as with the Binomial case in the next section), the model is being given three “opportunities” to apply Bayesian updating in estimating posteriors for each individual district, rather than a single one.
There is of the concern that the varying-intercepts model is being overfit in this case, but I will go forward with the analysis in this section nonetheless.
Fitting fixed-effects and varying-effects Bernoulli models (data from all three elections)
Having considered each individual election’s data and models, let’s fit a pair of models with one having fixed and the other having varying intercepts and slopes.
The fixed-effects model attempts to estimate a common intercept and common slope coefficients using observations across all cases - this is a “complete pooling” model where individual district-level elections are treated as all sharing a common underlying distribution for each parameter.
The verying-effects model attempts to estimate parameter values in a setting where there is some population-wide underlying set of values (as in the fixed-effects case), but where individual districts vary by differing degrees from these population-wide values, and estimation from each individual observation contributes to adjusted estimates for all other observations. This is the “partial pooling” model.
After evaluating the model fits and in-sample predictive performance, we can move on to sensitivity analysis (how do changes on the priors affect posterior distribution estimates?) and counterfactual analysis (how do changes in a given predictor’s values relate to changes in outcome prediction?).
Fitting the models
m1.allFull_fIfS <-
brm(generalWinD ~ 1 + incumbentD + incumbentR + primaryUnopD + primaryUnopR +
cmpgnMaxDisbursDvsR + primaryTotDvsR + log_popDensPerSqMi
+ medianAgeYr + pctEduGradProAge25plus +
pctInLaborForceUnemp + perCapitaIncome + pctWhiteAlone,
data = houseElectionDat_std, family = bernoulli(),
prior = c(set_prior("normal(0, 3)", class = "Intercept"),
set_prior("normal(0, 2)", class = "b")),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99),
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.allFull_fIfS")
m1.allFull_vIvS <-
brm(generalWinD ~ 1 + incumbentD + incumbentR + primaryUnopD + primaryUnopR +
cmpgnMaxDisbursDvsR + primaryTotDvsR + log_popDensPerSqMi
+ medianAgeYr + pctEduGradProAge25plus +
pctInLaborForceUnemp + perCapitaIncome + pctWhiteAlone +
(1 + incumbentD + incumbentR + primaryUnopD + primaryUnopR +
cmpgnMaxDisbursDvsR + primaryTotDvsR + log_popDensPerSqMi
+ medianAgeYr + pctEduGradProAge25plus +
pctInLaborForceUnemp + perCapitaIncome + pctWhiteAlone |
district2),
data = houseElectionDat_std, family = bernoulli(),
prior = c(set_prior("normal(0, 3)", class = "Intercept"),
set_prior("normal(0, 2)", class = "b"),
set_prior("cauchy(0, 2)", class = "sd"),
set_prior("lkj(2)", class = "cor")),
iter = 7000, warmup = 3000, chains = 3,
control = list(adapt_delta = 0.99),
thin = 2,
file = "C:/Users/Duane/Documents/Academic/Villanova/4. Fall 18/ProjectData/ModFits/m1.allFull_vIvS")Relative weights for fixed- vs. varying-effects models
In order to check how the two “full” models are estimate to perform in out-of-sample inference, let’s compare model weights using WAIC.
(m1.allFull_modWts <- model_weights(m1.allFull_fIfS, m1.allFull_vIvS,
weights = "waic") )## m1.allFull_fIfS m1.allFull_vIvS
## 8.998e-37 1.000e+00By a very wide margin, the varying-effects model is deemed the superior model for minimizing pointwise out-of-sample deviance - in other words, the varying-effects model is estimated to be less “wrong” for a given prediction on out-of-sample data. While I’ll be pursuing the most in-depth analysis on the varying-effects model, let’s continue to consider both models for now.
In-depth fit diagnostics
Just to confirm the models have fit well, I checked all density plots and trace plots for the various population-level and group-level effects - the fits are good (good stationarity around a concentrated set of estimated values, and good mixing of MCMC chains across the range of values).
Here are the model summaries (note that the R-hat values are all roughly 1.00, indicating good convergence):
summary(m1.allFull_fIfS)## Family: bernoulli
## Links: mu = logit
## Formula: generalWinD ~ 1 + incumbentD + incumbentR + primaryUnopD + primaryUnopR + cmpgnMaxDisbursDvsR + primaryTotDvsR + log_popDensPerSqMi + medianAgeYr + pctEduGradProAge25plus + pctInLaborForceUnemp + perCapitaIncome + pctWhiteAlone
## Data: houseElectionDat_std (Number of observations: 1305)
## Samples: 3 chains, each with iter = 7000; warmup = 3000; thin = 1;
## total post-warmup samples = 12000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept -0.48 0.31 -1.08 0.10 12000 1.00
## incumbentD 2.61 0.41 1.81 3.43 10115 1.00
## incumbentR -1.47 0.39 -2.25 -0.70 10151 1.00
## primaryUnopD 0.70 0.47 -0.22 1.63 7956 1.00
## primaryUnopR -1.35 0.50 -2.34 -0.37 7819 1.00
## cmpgnMaxDisbursDvsR 1.13 0.24 0.69 1.63 12000 1.00
## primaryTotDvsR 2.24 0.31 1.66 2.88 9489 1.00
## log_popDensPerSqMi 0.41 0.23 -0.02 0.86 12000 1.00
## medianAgeYr -0.21 0.20 -0.61 0.18 12000 1.00
## pctEduGradProAge25plus 0.13 0.39 -0.64 0.89 8201 1.00
## pctInLaborForceUnemp 0.71 0.23 0.27 1.16 12000 1.00
## perCapitaIncome 0.48 0.43 -0.35 1.32 7302 1.00
## pctWhiteAlone -0.40 0.30 -1.00 0.19 9686 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Fixed-effect model effective sample sizes are very large (the maximum possible is 12,000), and as noted the R-hat convergence diagnostic is good, suggesting the posterior estimates are fairly stable.
As might be expected, the mean estimates and 95% credible intervals for parameter posterior distributions vary, with many predictor posteriors overlapping both negative and positive values (and therefore also containing 0) - suggesting the association with the outcome variable generalWinD is uncertain and very plausibly small to nonexistent. Other parameters, such as incumbentD and incumbentR, are very plausibly nonzero and indicate relatively sizable association with a given probability of a Democratic candidate being elected.
summary(m1.allFull_vIvS)## Family: bernoulli
## Links: mu = logit
## Formula: generalWinD ~ 1 + incumbentD + incumbentR + primaryUnopD + primaryUnopR + cmpgnMaxDisbursDvsR + primaryTotDvsR + log_popDensPerSqMi + medianAgeYr + pctEduGradProAge25plus + pctInLaborForceUnemp + perCapitaIncome + pctWhiteAlone + (1 + incumbentD + incumbentR + primaryUnopD + primaryUnopR + cmpgnMaxDisbursDvsR + primaryTotDvsR + log_popDensPerSqMi + medianAgeYr + pctEduGradProAge25plus + pctInLaborForceUnemp + perCapitaIncome + pctWhiteAlone | district2)
## Data: houseElectionDat_std (Number of observations: 1305)
## Samples: 3 chains, each with iter = 7000; warmup = 3000; thin = 2;
## total post-warmup samples = 6000
##
## Group-Level Effects:
## ~district2 (Number of levels: 435)
## Estimate Est.Error l-95% CI u-95% CI
## sd(Intercept) 0.72 0.57 0.03 2.10
## sd(incumbentD) 4.70 1.36 2.12 7.43
## sd(incumbentR) 2.40 1.02 0.36 4.44
## sd(primaryUnopD) 1.26 0.94 0.06 3.47
## sd(primaryUnopR) 0.76 0.60 0.03 2.28
## sd(cmpgnMaxDisbursDvsR) 3.55 0.86 2.05 5.42
## sd(primaryTotDvsR) 0.85 0.64 0.03 2.38
## sd(log_popDensPerSqMi) 0.64 0.50 0.02 1.89
## sd(medianAgeYr) 0.76 0.54 0.03 1.98
## sd(pctEduGradProAge25plus) 0.65 0.49 0.03 1.84
## sd(pctInLaborForceUnemp) 0.73 0.55 0.03 2.01
## sd(perCapitaIncome) 0.59 0.47 0.02 1.73
## sd(pctWhiteAlone) 0.77 0.58 0.03 2.17
## cor(Intercept,incumbentD) -0.13 0.26 -0.60 0.40
## cor(Intercept,incumbentR) 0.01 0.25 -0.47 0.49
## cor(incumbentD,incumbentR) -0.25 0.21 -0.63 0.19
## cor(Intercept,primaryUnopD) 0.01 0.25 -0.48 0.48
## cor(incumbentD,primaryUnopD) -0.07 0.25 -0.53 0.42
## cor(incumbentR,primaryUnopD) 0.05 0.25 -0.43 0.52
## cor(Intercept,primaryUnopR) -0.01 0.25 -0.49 0.47
## cor(incumbentD,primaryUnopR) -0.01 0.25 -0.48 0.49
## cor(incumbentR,primaryUnopR) -0.03 0.25 -0.51 0.45
## cor(primaryUnopD,primaryUnopR) -0.03 0.25 -0.51 0.47
## cor(Intercept,cmpgnMaxDisbursDvsR) -0.04 0.24 -0.51 0.43
## cor(incumbentD,cmpgnMaxDisbursDvsR) -0.12 0.21 -0.50 0.31
## cor(incumbentR,cmpgnMaxDisbursDvsR) -0.03 0.23 -0.48 0.41
## cor(primaryUnopD,cmpgnMaxDisbursDvsR) -0.05 0.25 -0.51 0.43
## cor(primaryUnopR,cmpgnMaxDisbursDvsR) -0.01 0.25 -0.48 0.48
## cor(Intercept,primaryTotDvsR) 0.01 0.25 -0.48 0.49
## cor(incumbentD,primaryTotDvsR) -0.03 0.25 -0.50 0.46
## cor(incumbentR,primaryTotDvsR) 0.03 0.24 -0.44 0.51
## cor(primaryUnopD,primaryTotDvsR) -0.00 0.25 -0.49 0.48
## cor(primaryUnopR,primaryTotDvsR) -0.01 0.25 -0.50 0.48
## cor(cmpgnMaxDisbursDvsR,primaryTotDvsR) -0.07 0.24 -0.52 0.42
## cor(Intercept,log_popDensPerSqMi) -0.00 0.25 -0.48 0.48
## cor(incumbentD,log_popDensPerSqMi) -0.02 0.25 -0.49 0.45
## cor(incumbentR,log_popDensPerSqMi) -0.01 0.24 -0.47 0.46
## cor(primaryUnopD,log_popDensPerSqMi) 0.01 0.25 -0.47 0.48
## cor(primaryUnopR,log_popDensPerSqMi) 0.01 0.25 -0.49 0.49
## cor(cmpgnMaxDisbursDvsR,log_popDensPerSqMi) -0.01 0.25 -0.48 0.47
## cor(primaryTotDvsR,log_popDensPerSqMi) -0.02 0.25 -0.50 0.45
## cor(Intercept,medianAgeYr) 0.02 0.25 -0.48 0.49
## cor(incumbentD,medianAgeYr) 0.08 0.25 -0.41 0.55
## cor(incumbentR,medianAgeYr) 0.06 0.25 -0.43 0.53
## cor(primaryUnopD,medianAgeYr) -0.00 0.25 -0.49 0.48
## cor(primaryUnopR,medianAgeYr) -0.01 0.25 -0.49 0.48
## cor(cmpgnMaxDisbursDvsR,medianAgeYr) 0.01 0.24 -0.45 0.47
## cor(primaryTotDvsR,medianAgeYr) 0.00 0.25 -0.48 0.48
## cor(log_popDensPerSqMi,medianAgeYr) -0.02 0.25 -0.50 0.47
## cor(Intercept,pctEduGradProAge25plus) -0.01 0.25 -0.49 0.47
## cor(incumbentD,pctEduGradProAge25plus) -0.03 0.25 -0.51 0.46
## cor(incumbentR,pctEduGradProAge25plus) -0.02 0.25 -0.50 0.47
## cor(primaryUnopD,pctEduGradProAge25plus) -0.00 0.24 -0.47 0.47
## cor(primaryUnopR,pctEduGradProAge25plus) 0.00 0.25 -0.48 0.48
## cor(cmpgnMaxDisbursDvsR,pctEduGradProAge25plus) 0.01 0.25 -0.47 0.50
## cor(primaryTotDvsR,pctEduGradProAge25plus) 0.00 0.25 -0.49 0.48
## cor(log_popDensPerSqMi,pctEduGradProAge25plus) -0.01 0.25 -0.50 0.48
## cor(medianAgeYr,pctEduGradProAge25plus) -0.03 0.25 -0.51 0.46
## cor(Intercept,pctInLaborForceUnemp) -0.01 0.25 -0.51 0.47
## cor(incumbentD,pctInLaborForceUnemp) 0.05 0.24 -0.43 0.51
## cor(incumbentR,pctInLaborForceUnemp) -0.06 0.25 -0.53 0.43
## cor(primaryUnopD,pctInLaborForceUnemp) -0.02 0.25 -0.51 0.47
## cor(primaryUnopR,pctInLaborForceUnemp) -0.00 0.25 -0.49 0.47
## cor(cmpgnMaxDisbursDvsR,pctInLaborForceUnemp) -0.02 0.25 -0.48 0.46
## cor(primaryTotDvsR,pctInLaborForceUnemp) -0.02 0.25 -0.50 0.47
## cor(log_popDensPerSqMi,pctInLaborForceUnemp) 0.00 0.25 -0.48 0.48
## cor(medianAgeYr,pctInLaborForceUnemp) -0.01 0.25 -0.48 0.47
## cor(pctEduGradProAge25plus,pctInLaborForceUnemp) 0.02 0.25 -0.46 0.50
## cor(Intercept,perCapitaIncome) 0.00 0.25 -0.48 0.49
## cor(incumbentD,perCapitaIncome) -0.01 0.25 -0.49 0.48
## cor(incumbentR,perCapitaIncome) 0.00 0.25 -0.49 0.50
## cor(primaryUnopD,perCapitaIncome) 0.01 0.25 -0.47 0.49
## cor(primaryUnopR,perCapitaIncome) 0.01 0.25 -0.47 0.47
## cor(cmpgnMaxDisbursDvsR,perCapitaIncome) 0.02 0.25 -0.46 0.49
## cor(primaryTotDvsR,perCapitaIncome) 0.00 0.25 -0.49 0.49
## cor(log_popDensPerSqMi,perCapitaIncome) -0.02 0.25 -0.50 0.45
## cor(medianAgeYr,perCapitaIncome) -0.02 0.25 -0.50 0.47
## cor(pctEduGradProAge25plus,perCapitaIncome) -0.03 0.25 -0.52 0.45
## cor(pctInLaborForceUnemp,perCapitaIncome) 0.01 0.25 -0.48 0.49
## cor(Intercept,pctWhiteAlone) 0.01 0.25 -0.49 0.49
## cor(incumbentD,pctWhiteAlone) 0.02 0.24 -0.45 0.49
## cor(incumbentR,pctWhiteAlone) 0.04 0.25 -0.45 0.52
## cor(primaryUnopD,pctWhiteAlone) 0.02 0.25 -0.46 0.49
## cor(primaryUnopR,pctWhiteAlone) 0.00 0.26 -0.49 0.50
## cor(cmpgnMaxDisbursDvsR,pctWhiteAlone) 0.03 0.25 -0.44 0.51
## cor(primaryTotDvsR,pctWhiteAlone) 0.01 0.25 -0.45 0.49
## cor(log_popDensPerSqMi,pctWhiteAlone) 0.01 0.25 -0.47 0.49
## cor(medianAgeYr,pctWhiteAlone) 0.01 0.25 -0.47 0.49
## cor(pctEduGradProAge25plus,pctWhiteAlone) -0.01 0.25 -0.48 0.46
## cor(pctInLaborForceUnemp,pctWhiteAlone) 0.00 0.25 -0.48 0.48
## cor(perCapitaIncome,pctWhiteAlone) -0.00 0.25 -0.48 0.48
## Eff.Sample Rhat
## sd(Intercept) 2559 1.00
## sd(incumbentD) 2348 1.00
## sd(incumbentR) 2111 1.00
## sd(primaryUnopD) 3837 1.00
## sd(primaryUnopR) 4964 1.00
## sd(cmpgnMaxDisbursDvsR) 4534 1.00
## sd(primaryTotDvsR) 4434 1.00
## sd(log_popDensPerSqMi) 3724 1.00
## sd(medianAgeYr) 3484 1.00
## sd(pctEduGradProAge25plus) 4142 1.00
## sd(pctInLaborForceUnemp) 3281 1.00
## sd(perCapitaIncome) 4601 1.00
## sd(pctWhiteAlone) 4047 1.00
## cor(Intercept,incumbentD) 2420 1.00
## cor(Intercept,incumbentR) 3246 1.00
## cor(incumbentD,incumbentR) 3717 1.00
## cor(Intercept,primaryUnopD) 4065 1.00
## cor(incumbentD,primaryUnopD) 4876 1.00
## cor(incumbentR,primaryUnopD) 4658 1.00
## cor(Intercept,primaryUnopR) 4613 1.00
## cor(incumbentD,primaryUnopR) 4835 1.00
## cor(incumbentR,primaryUnopR) 4767 1.00
## cor(primaryUnopD,primaryUnopR) 4676 1.00
## cor(Intercept,cmpgnMaxDisbursDvsR) 3873 1.00
## cor(incumbentD,cmpgnMaxDisbursDvsR) 4354 1.00
## cor(incumbentR,cmpgnMaxDisbursDvsR) 3855 1.00
## cor(primaryUnopD,cmpgnMaxDisbursDvsR) 4145 1.00
## cor(primaryUnopR,cmpgnMaxDisbursDvsR) 4152 1.00
## cor(Intercept,primaryTotDvsR) 4080 1.00
## cor(incumbentD,primaryTotDvsR) 4611 1.00
## cor(incumbentR,primaryTotDvsR) 4782 1.00
## cor(primaryUnopD,primaryTotDvsR) 4837 1.00
## cor(primaryUnopR,primaryTotDvsR) 4881 1.00
## cor(cmpgnMaxDisbursDvsR,primaryTotDvsR) 4864 1.00
## cor(Intercept,log_popDensPerSqMi) 4891 1.00
## cor(incumbentD,log_popDensPerSqMi) 4349 1.00
## cor(incumbentR,log_popDensPerSqMi) 4759 1.00
## cor(primaryUnopD,log_popDensPerSqMi) 4815 1.00
## cor(primaryUnopR,log_popDensPerSqMi) 4555 1.00
## cor(cmpgnMaxDisbursDvsR,log_popDensPerSqMi) 4861 1.00
## cor(primaryTotDvsR,log_popDensPerSqMi) 4539 1.00
## cor(Intercept,medianAgeYr) 5008 1.00
## cor(incumbentD,medianAgeYr) 4594 1.00
## cor(incumbentR,medianAgeYr) 4650 1.00
## cor(primaryUnopD,medianAgeYr) 5047 1.00
## cor(primaryUnopR,medianAgeYr) 4528 1.00
## cor(cmpgnMaxDisbursDvsR,medianAgeYr) 4428 1.00
## cor(primaryTotDvsR,medianAgeYr) 5064 1.00
## cor(log_popDensPerSqMi,medianAgeYr) 5222 1.00
## cor(Intercept,pctEduGradProAge25plus) 4788 1.00
## cor(incumbentD,pctEduGradProAge25plus) 4515 1.00
## cor(incumbentR,pctEduGradProAge25plus) 4461 1.00
## cor(primaryUnopD,pctEduGradProAge25plus) 4710 1.00
## cor(primaryUnopR,pctEduGradProAge25plus) 4064 1.00
## cor(cmpgnMaxDisbursDvsR,pctEduGradProAge25plus) 4831 1.00
## cor(primaryTotDvsR,pctEduGradProAge25plus) 4449 1.00
## cor(log_popDensPerSqMi,pctEduGradProAge25plus) 4779 1.00
## cor(medianAgeYr,pctEduGradProAge25plus) 4884 1.00
## cor(Intercept,pctInLaborForceUnemp) 5015 1.00
## cor(incumbentD,pctInLaborForceUnemp) 4460 1.00
## cor(incumbentR,pctInLaborForceUnemp) 4810 1.00
## cor(primaryUnopD,pctInLaborForceUnemp) 5063 1.00
## cor(primaryUnopR,pctInLaborForceUnemp) 5038 1.00
## cor(cmpgnMaxDisbursDvsR,pctInLaborForceUnemp) 4978 1.00
## cor(primaryTotDvsR,pctInLaborForceUnemp) 4748 1.00
## cor(log_popDensPerSqMi,pctInLaborForceUnemp) 5174 1.00
## cor(medianAgeYr,pctInLaborForceUnemp) 4690 1.00
## cor(pctEduGradProAge25plus,pctInLaborForceUnemp) 5042 1.00
## cor(Intercept,perCapitaIncome) 4862 1.00
## cor(incumbentD,perCapitaIncome) 5066 1.00
## cor(incumbentR,perCapitaIncome) 4620 1.00
## cor(primaryUnopD,perCapitaIncome) 4596 1.00
## cor(primaryUnopR,perCapitaIncome) 4961 1.00
## cor(cmpgnMaxDisbursDvsR,perCapitaIncome) 4904 1.00
## cor(primaryTotDvsR,perCapitaIncome) 4669 1.00
## cor(log_popDensPerSqMi,perCapitaIncome) 4654 1.00
## cor(medianAgeYr,perCapitaIncome) 5048 1.00
## cor(pctEduGradProAge25plus,perCapitaIncome) 4382 1.00
## cor(pctInLaborForceUnemp,perCapitaIncome) 5089 1.00
## cor(Intercept,pctWhiteAlone) 4207 1.00
## cor(incumbentD,pctWhiteAlone) 4980 1.00
## cor(incumbentR,pctWhiteAlone) 4827 1.00
## cor(primaryUnopD,pctWhiteAlone) 4762 1.00
## cor(primaryUnopR,pctWhiteAlone) 4554 1.00
## cor(cmpgnMaxDisbursDvsR,pctWhiteAlone) 4603 1.00
## cor(primaryTotDvsR,pctWhiteAlone) 4913 1.00
## cor(log_popDensPerSqMi,pctWhiteAlone) 4883 1.00
## cor(medianAgeYr,pctWhiteAlone) 4872 1.00
## cor(pctEduGradProAge25plus,pctWhiteAlone) 4784 1.00
## cor(pctInLaborForceUnemp,pctWhiteAlone) 4892 1.00
## cor(perCapitaIncome,pctWhiteAlone) 4675 1.00
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept -0.50 0.63 -1.74 0.71 5127 1.00
## incumbentD 4.27 1.17 2.05 6.66 3728 1.00
## incumbentR -2.57 1.02 -4.71 -0.72 4069 1.00
## primaryUnopD 0.86 0.97 -1.09 2.68 4976 1.00
## primaryUnopR -2.38 0.99 -4.32 -0.38 4613 1.00
## cmpgnMaxDisbursDvsR 5.66 1.07 3.68 7.87 4926 1.00
## primaryTotDvsR 3.98 0.75 2.59 5.54 4663 1.00
## log_popDensPerSqMi 1.28 0.57 0.23 2.48 4901 1.00
## medianAgeYr -0.38 0.48 -1.37 0.56 4335 1.00
## pctEduGradProAge25plus 0.26 0.88 -1.43 2.01 4606 1.00
## pctInLaborForceUnemp 1.58 0.57 0.53 2.77 4873 1.00
## perCapitaIncome 0.66 0.91 -1.10 2.47 4551 1.00
## pctWhiteAlone -1.09 0.72 -2.54 0.28 5021 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).While the varying-effects model also has good convergence diagnostics, there are at least two major differences compared to the fixed-effects model:
There are many more parameters being estimated - this is because the “partial pooling” model is considering the divergence of each district-level value from the corresponding population-wide value for each parameter’s posterior.
Effective sample sizes are much smaller - this this is because I set the argument
thin = 2when fitting the model - this divides the number of post-warmup MCMC iterations to half of the “full” number, and setting “thin” > 1 reduces the computation time and file size of the model object. I did this because the complexity of the model (lots of observations and parameter estimations) resulted in a very large file (446MB, versus 891MB for the “default / thin = 1” equivalent).
In other words, the “true” model actually has approximately twice the effective sample size for each parameter, but our posterior estimates shown above are accurate to the “full” version of the model - I generated the “full” version of the model to confirm this, and just saved the “thinned” version. Estimates for the “full” version of the model had mean estimates and 95% credible interval bounds which differed from the “thinned” equivalent by less than 0.1 in all cases and for each parameter - I suspect this is simply due to my not having set a “seed” value for the MCMC random number generator.
Visualizing parameter posterior estimates
Comparing population-level effect estimates
First, let’s compare the population-level parameter posterior estimates between the fixed- and varying-effects models.
The varying-effects model has wider 90% credible intervals (indicating greater uncertainty in the posterior estimates) for each population-level parameter, and it does not have point estimates shrinking towards zero as I had expected. Among parameters with fixed-effects estimates farther from zero, the corresponding varying-effects estimates are even farther from zero - though the two 90% credible intervals largely overlap.
The two notable exceptions to this last note concerns cmpgnMaxDisbursDvsR (difference in maximum campaign-reported spending) and primaryTotDvsR (difference in total number of primary votes) - here, the varying-effects model is gives much more plausibility to a more positive effect.
The following plot shows the posterior slope estimates for the varying-effects model, which will be used to consider trimming the model.
The plot above suggests there are several parameters with 50% credible intervals overlapping zero (thus highlighting the plausibility such a parameter has a coefficient of zero, and therefore no marginal predictive value when the other predictors are accounted for). The most immediately plausible of these is the intercept, which I consider an essential component of the regression model. pctEduGradProAge25plus is the next likeliest candidate.
For the sake of convenience, I’ll evaluate the predictive merit of removing pctEduGradProAge25plus by comparing the WAIC weight of the “full” model versus that of the “reduced” model.
m1.allRed_vIvS <-
update(m1.allFull_vIvS, . ~ . -pctEduGradProAge25plus)
m1.all_modWts <- model_weights(m1.allFull_vIvS, m1.allRed_vIvS, weights = "waic")m1.all_modWts## m1.allFull_vIvS m1.allRed_vIvS
## 0.4318 0.5682Since the “full” and “reduced” models are allocated almost equal shares of the relative model weight, and the reduced model has only one fewer variable than the “full” model, I’ll just stick with the “full” model for now.
Reviewing group-level effects
I’m less sure of how to interpret parameters associated with varying intercepts/slopes, but I believe the above indicates that there is relatively little variability in slope coefficient values at the per-district level - that is, the estimated slope coefficient for a given district’s b_Intercept or b_perCapitaIncome value, say, is not likely to change very much from one year (in the 2012/’14/’16 data set) to another. It appears that this is not necessarily the case for the predictors incumbentD and cmpgnMaxDisbursDvsR, however.
The next plot’s a bit of an eyesore, but it seems useful to consider estimated correlations between the slopes of the various predictors.
It seems strongly plausible that the correlation between slopes for essentially every predictor variable pairing (including with the intercept) at the district level is 0. The estimated correlation between incumbentD and incumbentR has a point (mean) estimate markedly more different from zero than the others - though this particular interaction isn’t very meaningful since only a very small number of 2012 districts had both D and R incumbents due to post-2010 redistricting. I will consider (eventually) including interactions as additional predictors in the model.
These very-plausibly-zero correlation estimates may be a result of the brms package’s default to using non-centered parametrization, but I’m not sure.
Evaluating all-years Model 1 fit
Let’s assess how well m1.allFull_vIvS predicts the winner for in-sample data.
# prediction of each district (columns) for each post-warmup iteration (rows)
postPred_m1.allFull_vIvS <- posterior_predict(m1.allFull_vIvS)
# mean for each per-year district (column of postPred_m1.allFull_vIvS)
meanPred_m1.allFull_vIvS <- apply(postPred_m1.allFull_vIvS, 2, mean)
predVsObs.1v <-
tibble(district2 = houseElectionDat_std$district2,
year = houseElectionDat_std$year,
generalWinD = houseElectionDat_std$generalWinD,
meanPredWinD = meanPred_m1.allFull_vIvS )
nBins <- 100
cols <- c("darkred", "red",
"lightgrey",
"blue", "darkblue")
colGradient <- colorRampPalette(cols)
cut.cols <- colGradient(74)
cuts <- cut(predVsObs.1v$meanPredWinD, nBins)
names(cuts) <- sapply(cuts,function(t)
cut.cols[which(as.character(t) == levels(cuts))])
predVsObs.1v %>%
ggplot(aes(x = meanPredWinD, fill = cut(meanPredWinD, nBins))) +
geom_histogram(bins = nBins) +
facet_wrap(. ~ generalWinD) +
labs(title = "Posterior prediction check for m1.allFull_vIvS",
subtitle = "Faceted by Democratic candidate elected (0 / 1)",
x = "Mean proportion of model predictions, D candidate elected",
y = "Count") +
scale_fill_manual(values = cut.cols,
guide = F) +
theme_ds1()
The plot above, to my relatively inexperienced model-fitting eyes, seems excellent. The model doesn’t appear to have nearly as many “bad” misses where the mean of model-predicted outcomes infer the incorrect party’s success 75% or more of the time (relative to the single-year models - see Appendix 1). At the same time, the model seems to assign a fairly high mean predicted-outcome likelihood in a vast majority of cases - something on the order of 85%+ predicted-outcome mean proportions correctly inferring either the Republican or Democratic candidate’s winning.
For the sake of analysis, let’s create a quick “hit or miss” table similar to what was constructed for the individual year-fit models, but now defining a “miss” as a mean proportion of the election outcome assigning 40% or more of simulations to the incorrect party’s success.
# compute proportion of cases correctly assigned to generalWinD = 0 / 1
# using <= 40% & 0 // >= 60% & 1 as cutoff for acceptable accuracy
predVsObs.1v <-
predVsObs.1v %>%
mutate(hitOrMiss = if_else((generalWinD == 0 & meanPredWinD < 0.40) |
(generalWinD == 1 & meanPredWinD > 0.60), "hit",
"miss"),
missMargin = abs(generalWinD - meanPredWinD))
# data shown as counts
predVsObs.1v %>%
group_by(generalWinD, hitOrMiss, year) %>%
summarize(nObs = n())## # A tibble: 10 x 4
## # Groups: generalWinD, hitOrMiss [?]
## generalWinD hitOrMiss year nObs
## <int> <chr> <int> <int>
## 1 0 hit 2012 234
## 2 0 hit 2014 245
## 3 0 hit 2016 241
## 4 0 miss 2014 2
## 5 1 hit 2012 195
## 6 1 hit 2014 187
## 7 1 hit 2016 190
## 8 1 miss 2012 6
## 9 1 miss 2014 1
## 10 1 miss 2016 4# for which districts did the model "miss" most?
predVsObs.1v %>%
filter(hitOrMiss == "miss") %>%
arrange(desc(missMargin)) %>%
head(20) %>%
kable(align = "c") %>%
kable_styling(full_width = F, bootstrap_options = c("striped"))| district2 | year | generalWinD | meanPredWinD | hitOrMiss | missMargin |
|---|---|---|---|---|---|
| NE-2 | 2014 | 1 | 0.2558 | miss | 0.7442 |
| FL-7 | 2016 | 1 | 0.2662 | miss | 0.7338 |
| IL-12 | 2012 | 1 | 0.2773 | miss | 0.7227 |
| WV-3 | 2014 | 0 | 0.6642 | miss | 0.6642 |
| IL-17 | 2012 | 1 | 0.3787 | miss | 0.6213 |
| FL-13 | 2016 | 1 | 0.4290 | miss | 0.5710 |
| TX-23 | 2012 | 1 | 0.4695 | miss | 0.5305 |
| NY-24 | 2012 | 1 | 0.4817 | miss | 0.5183 |
| NV-3 | 2016 | 1 | 0.5110 | miss | 0.4890 |
| NY-18 | 2012 | 1 | 0.5190 | miss | 0.4810 |
| NH-1 | 2016 | 1 | 0.5602 | miss | 0.4398 |
| WV-3 | 2012 | 1 | 0.5673 | miss | 0.4327 |
| NV-4 | 2014 | 0 | 0.4010 | miss | 0.4010 |
Impressive. Across 1305 district-level elections (435 Congressional districts x 3 years of elections data), the model-predicted outcome ascribes a greater than 40% likelihood of the incorrect party’s success in 14 cases, with nearly half coming in 2012, the first year of 2010 Census-based redistricting. Interestingly, 11 of the 14 “misses” are for districts where Democratic candidates ultimately won - perhaps suggesting the model is in some way incorporating the fact that Republican candidates succeeded in a majority of districts each of the three years.
Considering whether the varying-effects model is overfit
I’m inclined to think the varying-intercepts model is appropriate for this analysis, but there is a definite possibility it is overfit to the data.
Let’s compare the full fixed-and varying-effects models once more, this time considering the predictions of the fixed-effects model and then the model density.
# prediction of each district (columns) for each post-warmup iteration (rows)
postPred_m1.allFull_fIfS <- posterior_predict(m1.allFull_fIfS)
# mean for each per-year district (column of postPred_m1.allFull_fIfS)
meanPred_m1.allFull_fIfS <- apply(postPred_m1.allFull_fIfS, 2, mean)
predVsObs.1f <-
tibble(district2 = houseElectionDat_std$district2,
year = houseElectionDat_std$year,
generalWinD = houseElectionDat_std$generalWinD,
meanPredWinD = meanPred_m1.allFull_vIvS )
nBins <- 100
cols <- c("darkred", "red",
"lightgrey",
"blue", "darkblue")
colGradient <- colorRampPalette(cols)
cut.cols <- colGradient(74)
cuts <- cut(predVsObs.1f$meanPredWinD, nBins)
names(cuts) <- sapply(cuts,function(t)
cut.cols[which(as.character(t) == levels(cuts))])
predVsObs.1f %>%
ggplot(aes(x = meanPredWinD, fill = cut(meanPredWinD, nBins))) +
geom_histogram(bins = nBins) +
facet_wrap(. ~ generalWinD) +
labs(title = "Posterior prediction check for m1.allFull_fIfS",
subtitle = "Faceted by Democratic candidate elected (0 / 1)",
x = "Mean proportion of model predictions, D candidate elected",
y = "Count") +
scale_fill_manual(values = cut.cols,
guide = F) +
theme_ds1()
# comparison of fixed- and varying-effects model predictions
ggplot() +
geom_density(data = predVsObs.1f,
aes(x = meanPredWinD), fill = color_set8[1], alpha = 0.4,
adjust = 0.2) +
geom_density(data = predVsObs.1v,
aes(x = meanPredWinD), fill = color_set8[2], alpha = 0.4,
adjust = 0.2) +
labs(title = "Comparing posterior prediction densities, 2016 models",
subtitle = "Fixed (dk blue) vs Varying effects (lt blue)",
x = "Mean proportion of model predictions, D candidate elected",
y = "Density") +
theme_ds1()
Since the density plot comparison shows an exact overlap between the two models, it seems as if the varying-effects model doesn’t contribute any additional inferential utility over the fixed-effects model.
I’ve gone through chapter 12 of Statistical Rethinking (which introduces and discusses partial pooling/multilevel models), and as best I can tell this is an appropriate application of the mechanism - so perhaps either model is actually just doing a good job of identifying that there are “Democratic”- and “Republican-friendly” districts.
Based on the definition of what constitutes a “miss” (described above), let’s compared how the two models fare.
# compute proportion of cases correctly assigned to generalWinD = 0 / 1
# using <= 40% & 0 // >= 60% & 1 as cutoff for acceptable accuracy
# now for the fixed-effects model
predVsObs.1f <-
predVsObs.1f %>%
mutate(hitOrMiss.f = if_else((generalWinD == 0 & meanPredWinD < 0.40) |
(generalWinD == 1 & meanPredWinD > 0.60), "hit",
"miss"),
missMargin.f = abs(generalWinD - meanPredWinD))
# adjust variable names for the varying-effects equivalent
predVsObs.1v <-
predVsObs.1v %>%
rename(hitOrMiss.v = hitOrMiss,
missMargin.v = missMargin)
# compare miss data shown as counts
predVsObs.1v %>%
bind_cols(predVsObs.1f) %>%
group_by(generalWinD, year) %>%
summarize(hit.f = sum(hitOrMiss.f == "hit"),
miss.f = sum(hitOrMiss.f == "miss"),
hit.v = sum(hitOrMiss.v == "hit"),
miss.v = sum(hitOrMiss.v == "miss") ) %>%
kable(align = "c") %>%
kable_styling(full_width = F, bootstrap_options = c("striped"))| generalWinD | year | hit.f | miss.f | hit.v | miss.v |
|---|---|---|---|---|---|
| 0 | 2012 | 234 | 0 | 234 | 0 |
| 0 | 2014 | 245 | 2 | 245 | 2 |
| 0 | 2016 | 241 | 0 | 241 | 0 |
| 1 | 2012 | 195 | 6 | 195 | 6 |
| 1 | 2014 | 187 | 1 | 187 | 1 |
| 1 | 2016 | 190 | 4 | 190 | 4 |
The two models appear to provide the same level of in-sample predictive capability, and I can’t come up with any better explanation than the fact that the two models are essentially different parametrizations of the same underlying model - in other words, there is no materially significant difference that makes the varying-effects model superior.
As a result, I’m relegating this section to Appendix 2 of this report and moving forward with the Binomial-outcome model section.
McElreath, Richard. Statistical Rethinking: a Bayesian Course with Examples in R and Stan. Chapman and Hall/CRC, 2016.↩
{https://cg-519a459a-0ea3-42c2-b7bc-fa1143481f74.s3-us-gov-west-1.amazonaws.com/bulk-downloads/index.html ; redirected from https://www.fec.gov/files/bulk-downloads/index.html }↩
https://factfinder.census.gov/faces/nav/jsf/pages/searchresults.xhtml?refresh=t↩
https://www.washingtonpost.com/blogs/the-fix/post/name-that-district-contest-winner-goofy-kicking-donald-duck/2011/12/29/gIQA2Fa2OP_blog.html?noredirect=on↩